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ABSTRACT:  The  motion  of  frontal  disturbances  In  the  atmosphere 
is  studied  based  on  several  nonlinear  models  proposed  by 
Stoker.   In  the  first  model,  the  air  Is  considered  to  be 
an  Incompressible  fluid  moving  over  a  plane  tangent  to  the 
rotating  earth.   The  fluid  consists  of  two  layers  and  the 
density  in  each  layer  is  assumed  to  be  constant.   The  hydro- 
static pressure  law  is  then  used  to  reduce  this  to  a  two  space 
dimensional  model.   The  boundary  between  these  layers  is  a 
contact  discontinuity  and  so  instabilities  may  occur  at  this 
frontal  surface. 

To  simplify  this  model, we  assume  that  the  dynamics  of 
the  perturbations  in  the  upper  warm  air  layer  can  be  neglected. 
In  this  case  only  the  motion  of  the  cold  air  need  be  studied. 
The  frontal  surface  intersects  the  horizontal  ground  in  a  curve, 
called  the  front,  which  is  a  free  boundary  for  the  cold  air. 
Following  the  procedure  of  Kasahara,  Isaacson  and  Stoker, 
we  make  a  numerical  study  of  this  model  using  generalizations 
of  the  Lax-Wendroff  scheme.   The  movement  of  the  front  is  based 
on  following  the  motion  of  material  particles  at  the  front. 
This  study  Indicates  the  development  of  the  occlusion  process 
from  an  initially  sinusoidal  frontal  pattern.   Thus,  we  show 
that  the  qualitative  features  of  the  occlusion  process  depend 
only  on  the  Coriolis  and  gravitational  forces  while  the 
thermodynamic  processes  can  be  ignored. 
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Various  initial  awu  u^undary  conditions  are  considered, 
and  a  comparl     f  their  effect  on  t  cess  is 

made.  In  all  cases,  the  cold  fr  faster  than  the 

warm  front  ,  -- -  -   relatively  sti^, .. 
exists  behind  i.  cycl 

pattern  al  t 

th;.  w,.^ — ....>,„  w;  ^^.v.^^  ^Lw....o  ^^sociated  wii;.  ■.  ..  ..^j, 

but  not  wl  ■ 

to  general  free  lems 

The  numerical  study  of  these  equations  is  still  • 
difficult  and  so  a  one  space  dimensional  model  is  int 

:-ical  comparison  with  the  two  dimensional  model  shows  that 

if led  theory  gives  many  of  the  important  character- 
frontal  motion  for  reasonable  lengths  of  time. 
The  one  dimensional  model  is  then  considered  for  a  seml- 
Inflr.  ■    "    ■         nstant  initial  and  boundary    ■  ■  •   ■. 

•St  order  hyperbolic  sys*  nded 

,_  -rles.  The  lowest  ordt :  ,  : ..._ :fy 

■ -' '  -UG  hyperbolic  system  of      '    . 

."Oducing  tl,        n 
invariants  .  by  Lax.   I. 

'.^onc    '•        '  ■  •       r.  ...       ^    '     "  ■  • 

are  given.  Wt.         continuou:  ,  they  can 

y  for 

>'•  .-•'- _..  .-        ,.  ,.  .^^   nonhi  ■ 

which  ca:  tly  for  the  terms  in 
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the  expansion.   Comparison  of  the  series,  through  second 
order  terms,  with  the  numerical  solution  of  the  original 
system  shows  close  agreement  away  from  the  boundary. 
These  techniques  can  be  used  for  all  nonhomogeneous 
quasi-linear  systems  where  the  solution  to  the  homogeneous 
system  is  known. 
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Introduction 

In  meteorology  It  Is  known  that  the  weather  In  the 
middle  lattltudes  Is  conditioned  to  a  considerable  degree 
by  events  associated  with  the  propagation  of  wave-like 
disturbances  on  a  discontinuity  surface  between  warm  and 
cold  air  masses  In  the  atmosphere.   The  Intersection  of  the 
discontinuity  surface  with  the  earth  Is  known  as  a  front. 
In  this  paper  we  will  be  Interested  In  following  the  motion 
of  these  fronts. 

The  Importance  of  nonlinear  effects  In  the  dynamical 
equations  of  the  cold  front  were  first  pointed  out  by  Freeman 
(1952)  [6],  Abdullah  (ig'^g)  [1],  and  Tepper  (1952)  [l8]  who 
applied  the  method  of  characteristics  to  solve  nonlinear 
equations  In  one  space  dimension.   Later  Stoker  (1953)  [l6] 
developed  a  two-layer  model  for  the  motion  of  a  frontal 
surface.   Whltham  (1953)  [20]  made  a  qualitative  study  of 
this  model  which  strongly  indicated  that  the  evolution  of 
frontal  cyclones  might  well  follow  the  pattern  that  leads 
to  the  occlusion  effect  observed  in  nature.   Stoker  together 
with  Kasahara  and  Isaacson  (196^)  [8,9]  made  some  numerical 
calculations  with  a  one-layer  model  and  found  the  beginnings 
of  the  occlusion  effect.   One  of  the  objects  of  this  paper 
is  to  continue  and  extend  the  investigations  of  Kasahara, 
Isaacson  and  Stoker  (abbreviated  as  K.I.S.). 

In  Chapter  II  we  will  present  several  finite  difference 
schemes  for  dealing  with  the  problem  under  various  initial  and 
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boundary  ----'•■  — .   --  -^~    -■--    t.-.-v,  theory  we  assume 
that  the  the  warm  air  is  not  affected  by  th- 

:  air  ;■  remains  in  its 

Initial  state.   This  *.  -  —.   the  followinr  --'---'--  -*' 

Dps  a  bulge  the  wai'        n 
compe:.  /a  slight  change  in  its  vertical 

-"'-■^-nent  ?.•-''  ■•''th  hardly  any  change  in  the  horizontal 

.-.ents  u  and  v.   Since  we  are  ignoring  changes  in  the 
vertical  velocity  in  any  case,  it  is  reasonable  to  assume 
that  the  warm  air  is  unaffected  by  these  movements  of  the 
front.   However,  in  the  cold  air  It  is  evident  that  relatively 
large  changes  in  the  components,  u  and  v,  of  the  velocity  may 
be  needed  to  move  around  such  a  frontal  disturbance. 

In  the  one  layer  model  the  discontinuity  surface 
between  the  warm  and  cold  air  masses  Is  a  free  surface.   The 
motion  of  this  free  surface  is  determined  solely  by  the 
dynajnics  of  the  cold  air.   However,  in  the  two  layer  model 
this  free  surface  becomes  a  contact  discontinuity  separating 
the  warm  and  cold  air  masses.   It  is  well  known  that  the 
occurrence  of  contact  discontinuities  leads  to  Instabilities 
both  physically  and  numerically  (e.g.  see  Richtmyer  [l'^]). 
In  Chapter  III  we  will  give  some  preliminary  results  using 
the  two  layer  model. 

In  Chapter  II  we  will  discuss  the  one  layer  model. 
We  are  able  to  continue  the  solution  for  longer  times  that 
K.I.S.  v;  illows  us  to  achieve  greater 
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insights  Into  the  formation  of  occluded  fronts.   It  Is 
planned  to  Incorporate  this  theory  Into  a  large  scale  model 
for  the  circulation  of  the  atmosphere.   The  motion  of  the 
fronts  could  then  be  followed  by  using  a  refined  mesh  in  the 
neighborhood  of  the  front.   Ciment  [4]  has  shown  that  such 
difference  schemes  with  uneven  mesh  spacings  are  stable, 
in  the  linear  case,  for  dissipative  schemes  as  the  Lax- 
Wendroff  method. 

Even  this  simplified  model  is  very  complicated  and 
the  equations  can  only  be  solved  numerically  and  with 
considerable  difficulty.   Therefore,  Stoker  [l6]  Introduced 
another  model  which  contains  only  one  space  dimension  but 
four  dependent  variables.   This  is  done  by  introducing  a 
long  wave  theory  in  the  horizontal  plane.   This  model  is 
derived  in  Chapter  1  since  the  hypotheses  of  this  model 
can  be  used  to  derive  boundary  conditions  for  the  two 
dimensional  model. 
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ch; 


a.  •!__ -  _v 

The  earth  i  rs 

of  air, 

to  th'-  ,—  "-.'.    '-'■'    --•-  ------.   The  density  of  th'   ■" - 

wl'  assumed  to  be  constant, 

Th:      ut  this       we  consider  a  tangent  plane  apf : 
tic-.  *- -    *-':e  earth.   The  x-axis  is  the  normal  to  this  plane, 
:tive  direction  being  directed  away  from  the  earth. 

.  and  y  axes  are  in  the  plane  of  the  earth  with  the 
positive  x  axis  directed  to  the  East  and  the  positive  y  axis 
to  the  North.   Thus,  we  are  ignoring  the  sphericity  of  the 
earth  even  though  we  shall  not  ignore  its  spin.   Because  of 
the  earth's  spin  we  shall  assume  that  the  coordinate  system 

tating  about  the  z  axis  with  a  constant  angular  velocity 
Q   =   iji   sin  0  =  ^  f,  where  w  is  the  angular  velocity  of  the 
'-ci:'.:.,  <\)    is  the  latitude  of  the  origin  of  the  system  and 
f  is  the  (constant)  Coriolis  parameter. 

Initially  th»=»  cold  air  lies  in  a  wedge  under  the  warm 
air  and  the  disc  m,  _  Li.uity  surface  between  the  two  layt^  ;  ..  is 
Incl'  'he  horizontal  as  shown  in  Figures 

nulty  surface  is  inclined  with  respect 
-  ,  »..  .  ..............  .........  .,..:,.,...,  of  the  Coriolis  ef  f'--'- . 

'er,  t'  quite  small,  of  the  order 

'>y  in  the  warm  air  is  that  of  th'     .ailing 
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westerlies  while  In  the  cold  air  the  velocity  Is  approxi- 
mately that  of  the  polar  wind.   Thus  across  the  discontinuity 
surface  there  are  tangential  discontinuities  in  the  velocities 
as  well  as  jumps  in  the  densities. 

Following  the  presentation  of  Stoker  [16]  we  begin 
with  the  Euler  equations  for  fluid  dynamics. 

/  N     du      8p  ,    „ 

^1-1)   ^^)   P  dt  =  -  9f  +  P  ^(y)     dt  =  9t  -^  ^  9Y  -^  ^  97  ^"  9^ 

/■     \  dw      9p  .    T-, 

^^^   P  dt  =  -  it  +  P  ^z)  • 

F  is  due  to  both  gravitational  and  Corlolis  forces. 
In  addition  we  assume  that  the  air  Is  incompressible  and 
of  constant  density  in  each  layer  and  so  we  have 

(d)     ^+   p.+   pl=    0    . 
9t    9y    9z 

We  shall  call  this  problem  I. 

In  dynamic  meteorology  a  standard  assumption  is  the 
hydrostatic  pressure  law.   This  states  that  the  vertical 
accelerations  of  the  Corlolis  term  in  the  third  equation 
of  (1)  can  be  Ignored  and  hence 

(1.2)  II   =   _p  g 

Thus,  on  purely  kinematic  grounds  we  have 

u  =  u(x,y,t)  ,    V  =  v(x,y,t),    w  =  0. 
Thus  we  have 
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'ic,-:i.z  1 


cold 


waroi 


EA3T 


h'  (  X I  y  ,  t ) 


V  =, 


h(K,y,t) 


w^x^yyjis^m^^TX^^^^?: 


NO'^TH 
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(a)             u,  +UU      +VU      =--P      +F/S 

t  X             y             p      X           (x) 

(1.3)       (b)             v^  +   ^v^   +    vvy    =   -   i-Py    +   F(y) 

(c)  u^    +      ^    ^    °    ' 
where 


F/    N    =    203   sin    ())    •    V   =    fv 
P/    >,    =-2u)    sin    d)    •    u   =-fu 

(y) 


f  a  constant 


Before  continuing  we  wish  to  Introduce  notation  to 
distinguish  between  the  warm  and  cold  air  layers.   Thus, 
from  now  on  primed  variables  e.g.  u',  v'   ref  r  to  the  warm 
air  mass  while  unprlmed  variables  e.g.  u,  v   refer  to  the 
cold  air. 

We  now  Integrate  equation  (1.2)  -^  =   -p    g.      Let 

a  Z 

h  =  h(x,y,t)  be  the  vertical  height  of  the  discontinuity 
surface  between  the  warm  and  cold  air  and  let  h'  =  h'(x,y,t) 
be  the  total  height  of  the  warm  air.   Assume  p'  =  0  at  the 
top  of  the  warm  air.   We  then  have 

p ' (x,y  ,z,t)  =  p  'g(h'-z) 
(1.^) 

p  (x,y,z,t)  =  p'g(h'-h)  +p  g(h-z)  , 

where  we  have  used  the  continuity  condition  p'  =  p  at  z  =  h. 
Substituting (1 . 4)  Into  (1.3)  we  get 

(a)  u^  +  uu^  +  vuy  =  -g[^  h-   +  (1-  ^)h^]  +  fv 

(b)  v^  +  uv^  +  vvy  =  -glf   hy   +  (1-  ^)hy]  -  fu 


(1.5) 


(c)      u   +  V   =  0 
X     y 
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(d)  u^  +  u'u^  +  v'u   =  -  gh^  +  fv' 

til        I 

(e)  V   +  u'v   +  v'v   =  -  gh   -  fu' 

t       A       y        y 

(f)  u'  +  V.' 

The  kinematic  at  the  free  surfaces  z  =  h 

,  ;  .'  =  h  are  f^^^"^^  ^  °'  HT^^"^'^  ^  °-  2^"°^  If  "  ° 
because  we  are  ignoring  vertical  velocities,  we  can  rewrite 

these  equations  using  partial  derivatives  instead  of  particle 

derivatives  and  get 

uh   +  vh   +  h.  =  0 
X     y    t 

u'h   +v'h   +  h.  =  0 
X     y     t 

u'h'  +v'h'  +  h^  =  0  . 
X     y    t 

Comoining  these  equations  together  with  (1.5c,f)  we  get 
h^  +  (hu)„  +  (hv)   =  0 


(1.6) 


t   '   'x      y 
(h'-h)^  +  [(h'-h)u]^  +  [(h'-h)v]y  =  0 


Physically  equations  (1.6)  are  simply  the  equations  for  the 
conservation  of  mass  in  both  the  warm  and  cold  air.   If  we 
substitute  (1.6)  into  (1.5)  in  place  of  equations  (1.5c,f) 
we  arrive  at  problem  II. 


(a)     u^  +  uu^  +  vuy  =  -g[^  h^  +  (1-  e-)h^]  +  fv 


P_  V.        X   M  _  £_ 
P 

II  p» 


(b)      V^  +  UV^  +  VVy  =  -g[^  hy  +  (1-  ^)hy]  -  fu 
(1.7), II 


'  •)      h^  +  (hu)^  +  (hv)y  =  0 


»  t  t  T 

(d)  u,  +  u'u   +  v'u   =  -gh   +  fv' 

t      X      y    ^  X 

f       t       t       I 

(e)  V,  +  u'v   +  v'v   =  -gh   -  fu' 

t     X     y    ^  y 

(f)  (h'-h)^  +  [(h'-h)u]^  +  [(h'-h)v]y  =  0 

System  II  has  an  exact  solution  corresponding  to  an 
Initial  state  of  a  stationary  front 

(a)  u'  =  u' 

v'  =  0 

h'  =  -  i|-  (y+K-^) 

(1.8) 

(b)  u  =  u 


V   =  0 

g(l-  ^)   P  ^ 


b .   Two-dimensional  -  Single  layer  theory 

The  two  layer  theory  just  given  Involves  several 
numerical  difficulties.   First,  there  are  six  dependent 
variables.   Since  the  problem  Involves  two  space  dimensions 
we  require  large  matrices  to  record  all  the  known  values 
of  these  variables  at  all  the  mesh  points  and  for  some  time 
step.   Because  of  this  the  mesh  must  be  fairly  coarse. 
Another  difficulty  is  the  high  sound  speed  in  the  warm  air. 
The  largest  sound  speed  of  the  two  layer  model  is  of  the 
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ord-     yS^- —   since  — '^'  1  and  h'  >>  h  this  sound  speed 
'   P  P  - 

is  considerably  larger  than  that  of  the  cold  air  which  Is 


of  *■•  -   rder  of  /g(l-  ^)h  •   By  the  Courant-Freidrichs-Lewy 
theory  [5]  the  maximum  time  step  allowable  for  stability  is 
inversely  proportional  to  the  maximum  sound  speed.   Thus  the 
time  s*--"-  must  be  very  small  and  so  the  solution  would 
require  large  amounts  of  computer  time.   However,  since  most 
of  the  dynamics  is  in  the  cold  air  the  sound  speed  of  physical 
relevance  is  that  of  the  cold  air.   Thus  the  small  time  step 
is  necessary  for  numerical  rather  than  physical  reasons. 
A  final  and  more  serious  difficulty  is  that  the  discontinuity 
surface  between  the  warm  and  cold  air  masses  Is  a  contact 
discontinuity.   Thus  Rayleigh  and  Helmholtz  Instabilities 
may  occur.   Even  in  the  absence  of  physical  instabilities 
various  numerical  Instabilities  occur  in  the  neighborhood  of 
contact  discontinuities. 

For  these  various  reasons  it  becomes  desirable  to 
eliminate  the  influence  of  the  warm  air  on  the  cold  air. 
The  physical  Justification  of  this  assumption  was  discussed 
in  the  Introduction.   Problem  III  is  gotten  by  assuming 
that  the  perturbations  in  the  cold  air  do  not  affect  the 

he  initial  conditions  (1.8a)  hold  for  all 
time.   Si.     uting  (1.8a)  into  system  II  we  have  system  III. 

^  +  ""x  -^  ^%  +  s(i-  r")^  =  ^^ 

(1.  ^         V.  +  uv^  +  vv,,  +  g(l-  ^)h„  =  f(u-  ^  d') 

u      A      y  p    y  P 

h^.  +  (hu)^  +  (hv)   =  0  . 
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In  this  system  there  are  three  dependent  variables ,  the 

2        o ' 
sound  speeds  are  c   =  g(l-  ■^^^— )h  and  the  front  is  a  free 

surface.   It  was  this  problem  that  Kasahara,  Isaacson 

and  Stoker  discussed  in  their  paper. 


c .   One  dimensional  model 

Stoker  [16]  developes  problem  IV  as  an  approximation 
based  on  assuming  waves  of  wave  lengths  that  are  large  compared 
to  the  typical  north-south  lengths.   Let  n(Xjt)  be  the 
displacement  of  the  front  from  the  rigid  wall   y  =  6  as  shown 
in  Figure  2.   The  velocity  v(x,y,t)  is  assumed  to  be  zero 
at  the  rigid  wall  and  then  increase   linearly  to  its  value 
at  the  front  y  =  6-ri  (x,t)  where  it  is  denoted  by  v(x,t). 
h  is  zero  at  the  front  and  increases  linearly  in  y  till  the 
rigid  wall  y  =  6.   Also,  we  denote  the  intersection  of  the 
discontinuity  surface  z  =  h(x,y,t)  with  the  plane  y  =  6  as 
h(x,t).   Finally  we  assume  that  u  is  independent  of  y  and 
hence  u  =  u(Xjt).   Thus  we  have 

(a)  h(x,y,y)  =  lzfi±^ll|^h(x,t) 
(1.10) 

(b)  v(x,y,t)  =  ^  (^~  ^^    v(x,t) 


In  addition  we  assume  that  a  particle  that  starts  on  the  front 
y  =  6-ri  (x,t)  remains  on  the  front  and  so 

(c)     v(x,t)  =  -  (n^  +  un^)  . 


This  is  slightly  different  from  Stoker's  definition  and 
corresponds  to  6- n  in  his  book. 
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FIGURE  ? 


yrr 


n(x,t)  Is  the  distance  from  the  front  to  the  wall  y  =  6 
h(x*,t)  is  the  height,  h(x,y,t)  at  the  wall  y  =  6. 
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We  now  integrate  the  equations  in  system  III  from  y  =  6-  ri 
to  y  =  6  and  obtain  system  IV.   To  show  how  this  is  done 
we  calculate  several  integrals  explicitly. 


6- 

-  V                         6- 

-  n 

c 

• 

S                                c 
[      V   ay   =   V 

S 

[y-(6-  n  )]  dy  =  -^  h  n 


(6-y)  dy   =  -^  y  r\ 
6- n         6- n 

We  now  differentiate  these  formulas  keeping  in  mind  that 

the  limits  of  integration  involve  n  which  is  a  function 

of  X  and  t .   Thus  we  have 


So 


3x 


V  dy  = 


v^  dy  -  v(x,6-n  ,t)  3^(5-11  ) 


6-n 


6-n 
6 


V   dy  +  v  n 

X    "^         X 


6 
r 


6-Tl 

6 


^x  ^y  =  97 


V  dy  -  vn^  =  11^^  (vn  )  -  vn^ 


6-  n 


6-  n 


or 


I  v^  n  +  I  V  n  ^ 


-  V  n 


V   dy  = 
X   "^ 


I  ^x^ 


I  ^^x 


6-n 


Similarly 
5 


h^  dy  =   I  K^  n  +  I  h  n  ^ 


6-n 


since  h(x,6-n  jt)  =  0, 
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I   ^  ^y  =   I  ^  ^  -  I  ^  ^  t 


6-n 
5 


I  ^t  ^y  =  I  i^t  ^  "■  I  ^  ^  t 


6-  r 


Also, 


j   (hv)y  dy  = 
6-n 

6 

I  (hu)^  dy  = 
6-  n 


hv 


6-n 


=  0   since  h(x  ,6- p  ,t )=v (x ,6 ,t )  =  0 


k^- 


h   dy] 


6-n 


I  (V  +  nu^)n  +  |hun^ 


=   2  ^^^^x  ^  '*'  2  ^^  ^  X  * 


>w  integrate  system  III  with  respect  to  y  from 

6-n  to  6  and  then  divide  the  resulting  equations  by  n  . 
k  =  g(l  -  ^^—)    and  we  have 

"^t  -^ "%  ^  I  ^^  -^  1 1^  ^  X  =  I  ^^ 

(1.11)         v^  .uv^  -  Vn-^  -%^  =  T-^-  ''^^-   f^'^ 


If  we  als-  '--^  ;)-  .-.quation  (1.10c)  and  use  it  to  simplify 
the  ab      mations  we  arrive  at  system  IV. 
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(a)  u^  +  uu^  +  I  kh^  +  I  ^  n^  =  I  fv 


(1.12), IV 


(b)   h^  +  (hu)^ 


(c)   V,  +  uv 

t        X 


hv 

n 


^^  -2f(u-  ^  a-) 

n        p 


(d)  '^t  "*"  ^  '^x 


=  -  V 


It  is  convenient  to  change  variables  by  introducing  the 

2    1 
sound  speed  c   =  p-  kh ,   we  then  have 


(a)   u,  +  uu  +  2cc   +  —  n 
'    t      X      X    n    X 


Ifv 


(1.13) 


(b )  c,    +   ^  u     +  uc 
t    2   X     X 


(c)   V,  +  uv 
t      X 


1_  cv_ 
2   n 


_^  _c_  _  2f(u_  Pj.  u,) 


(d)   '^t  "*"  '^  '^x 


=  -  V 


The  scheme  is  discussed  in  detail  by  Turkel  [191. 
Numerical  solution  of  this  model  shows  qualitative  agreement 
with  the  single  layer  theory.  Problem  III,  through  the  first 
eight  hours.   This  system  is  also  considered  in  a  semi- 
infinite  domain  and  a  formal  perturbation  series  solution 
is  obtained.   The  lowest  order  term  is  analyzed  by  use  of 
the  Lax  theory  of  Rlemann  Invariants  for  systems  of 
equations  [12].   The  first  few  terms  of  this  series  shows 
close  agreement  with  the  numerical  solution  of  the  system 
of  equations. 
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CHAi  --..  --.   ....  NAL  -  SINGLE  LAYER 


I:.-  ._,-. 


In  this  cha;  juss  problem  III, 

•  ^.  -^  ions  (1.9).   -     .3del  Is  a  two-space  dlme;.___; 

'em  in  which  we  assume  that  the  warm  air  remains  in  its 
initial  state  for  all  time  and  hence  we  can  concern  ourselves 
with  the  motion  of  the  cold  air  only.   We  assume  that  the 
cold  air  lies  above  a  region  D  of  the  x-y  plane  and  that  this 
region  is  bounded  on  three  sides  by  straight  lines  and  on 
the  fourth  side  by  a  curve  C(t)  as  illustrated  in  Figure  3. 
The  warm  air  lies  above  the  cold  air  in  three  dimensional 
space  and  occupies  the  entire  rectangle  R  in  the  x-y  plane. 
The  curve  C  is  the  line  where  the  cold  air  ends  i.e.  where 
h  =  0.   According  to  the  Glossary  of  Meteorology  (American 
Meteorological  Society  1959)  the  intersection  of  the 
discontinuity  surface  between  the  cold  air  and  warm  air 
with  the  earth's  surface  is  called  a  surface  front.   However 
we  shall  follow  popular  usage  and  call  C  the  front. 

In  this  discussion  we  confine  ourselves  to  a  region  D 
callei  "'  '    ■^  '  air  mass  domain.   The  curve  C(t):(x  "',   't)) 
is  which  h  =  0.   The  i 

velocity  (Up,Vp)  of  the  parti 


.  t-  ,.   r» 


ront,  ani  '      '"   following  conditions  hold 
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h  =  0 

where  d/dt  is  the  particle  derivative. 

We  must  still  prescribe  u,v  and  h  along  the  other 
three  boundaries.   For  simplicity  we  shall  assume  that  u,v 
and  h  are  periodic  in  the  space  variable  x  with  a  period 
equal  to  the  distance  between  the  east  and  west  boundaries. 
Since  the  atmospheric  motion  on  the  earth  is  periodic  in 
the  longitude  this  condition  is  correct  on  a  global  scale. 
Furthermore,  if  occlusion  takes  place  in  a  relatively  small 
region  centered  in  the  period,  we  may  by  varying  the  size 
of  the  period  determine  the  influence  of  the  periodicity 
condition  and  calculations  are  made  with  this  object  in  view. 
We  will  also  consider  one  case  in  which  the  boundaries  are 
so  far  apart  as  to  have  no  influence  on  the  occlusion  process 
and  so  the  domain  can  be  considered  as  Infinite  in  the  x 
direction. 

Along  the  northern  boundary  we  assume  that  the  normal 
component  of  the  velocity,  i.e.  the  y  component  of  velocity, 
v,  vanishes  for  all  time.   This  is  in  agreement  with  the 
assumption  made  in  deriving  the  one   dimensional  model, 
problem  IV,  so  that  a  meaningful  comparison  between  the 
two  models  can  be  made. 
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A  complete  formulation  of  the  problem  requires 
appropriate  initial  data  at  the  time  t  =  0  for  u,v,h  in 
the  cold  air  mass  domain  D.   It  is  also  necessary  to  specify 
the  initial  shape  of  the  curve  C  and  the  values  u,v  along 
the  front.   The  curve  C  is  assumed  to  be  sinusoidal 
initially.   The  thickness  h  of  the  cold  air  is  assumed  to 
vary  linearly  in  y  and  to  be  equal  to  zero  along  the  front; 
thus  h  varies  sinusoidally  in  x .  u  and  v  are  assumed  to 
be  those  of  a  steady  state  solution  of  equation  (1.9). 
Since  the  essential  change  made  from  the  steady  state 
solution  occurs  only  in  the  shape  of  the  front  our  initial 
conditions  can  be  considered  a  perturbation  of  the  steady 
state  solution  of  equation  (1.9). 

Two  different  sets  of  initial  conditions  are 
taken,  as  follows. 

Case  1 

^C  ^  *"2  ~  ^1  ^°^  ^-^T"  ^^ 
u   =   u 

V   =   0 

h  =   (y-y(.)H  ,  y^  <  y  <  Y 

where        H  = —   (—  u'  -  u)  ,   -2-^  u '  >  u  . 

g(l-  §-)  ^  " 

X,  Y  denote  the  lengths  of  the  sides  of  the  rectangle  R  in 
the  X  and  y  directions  respectively  and  u  is  the  x-velocity 
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c:  ...J   warm  air.   As  a  second  case  we  assume  a  physically 

more  relevant  condition,  i.e.  that  initially  the  wind 

Is  geostrophic.   This  condition  means  that  initially  there 

is  no  acceleration  in  either  the  x  or  y  terms.   Thus  we 

choose  our  initial  velocities  u  and  v  so  that  -r^  =   0   and 

at 

^  =  0 
dt   "• 

Case  2 

y  ~y 
h  =   (rr— ^)(Y-b)H 

g(l-  ^)    .,       . 

f     9y   p 

^  f     3x 

where   b  =  C-,  +  Cp  ,    H,Y  and  y   as  in  Case  1; 
with  these  initial  conditions  we  have  at  t  =  0 


§=  u^  +  uu^  +  VUy  =  -g(l-  f)h^   +  fv  =  0 

fj  =  v^  +  uv^  -.  vv^  =  -g(l-  f)h^   -    f(u-  fu^)    =    0. 


Since  we  have  as  an  initial  condition  u  =  u  the  bulge 
in  the  front  will  begin  to  move  to  the  right  and  away  from 
the  center  of  the  region.   We  wish  to  keep  this  bulge  as 
centered  as  possible,  so  that  the  occlusion  pattern  will  not 
be  affected  by  the  periodicity  condition.   It  is  thus 
convenient  to  introduce  a  moving  coordinate  system  so  that 
initially  the  front  is  stationary  in  the  x  direction.   Since 
*  rlodicity  condition  will  be  imposed  in  this  moving 
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coordinate  system  it  is  possible  to  thus  minimize  the 
effects  of  the  periodicity  condition  on  the  occlusion  process 
by  forcing  the  bulge  in  the  front  to  remain  near  the  center 
of  the  domain.   At  the  same  time  we  also  introduce  dimension- 
less   variables.   We  therefore  introduce  the  following  new 
variables .  ■• 

_   t  ,  _  At 

'^  "  At  '  ^  "  A? 

^  As  ^   As 

(2.1)      .        _ 

u  =  A.(u-u)  V  =X  V 

h  =  A^gd-  fi^)h 

where  u  is  the  constant  initial  velocity,  in  case  I,  of 

the  cold  air  mass.   At,  As   denote  units  for  time  and  length 

respectively.   We  also  introduce  the  following  parameters: 

F  =  f  At  G  =  FA(^  u'  -  u)  . 

P 

With  these  new  variables  the  equations  become 


/\    /V  >N    ^ 


(a)  u     +   uu^   +   vu      +   h^   =   Pv 

(2.2)       (b)  v^    +    uv^    +VV      +h      =-Fu+G 

(c)  h      +    h(u_+    V       )    +    uh>    +    vh      =    0 


with  initial  conditions 
Case  I 

(d)      u(0,C,  r,  )  =  0 
v(0,^,  n  )  =  0 


^2    ^1 


h(0,C,  n  )  =  G(  n  -  n  ^)  where  "^  q=  -^  -  Xs 


,2t\    As  ^-^ 
'COS( — j^ —  ?) 
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Case  II 


(d)'    h(0,C,n)=G(n-n  (,)  (N-b)/(N- n  ^) 
u(0,^,n  )  =  ^  (G  -  14.(0, C,n  )) 


and 


F 
1  3h 


3  n' 


v(0,C,Ti  )  =  ^  If  (0,C,Ti  ) 
Y 


where   N  = 


As 


b  = 


As 


3h 
3  n 

p/N    -    bx 

3h 

G(N-b)(N-n)      ^"^c 

3? 

(N-.,)2            n 

9^C 

2TfC,                   „       . 

.1        r.in     f2Tr_As     p 

3C 


X 


as  before. 


Notice  that  v(0,C,N)  =  0  and  so  the  initial  conditions 
match  the  boundary  conditions  at  the  northern  boundary  as 
well  as  at  the  other  boundaries.   The  boundary  conditions 
at  the  northern,  eastern  and  western  boundaries  are 

(e)     v(t,^,N)  =  0 

periodicity  in  the   C  direction. 

While  the  boundary  conditions  along  the  front  are 


(f) 


h(T,C,  n  ^)  =  0 

dXp 

dT    ^C 


dT 


=  yVp  -  Vh  +  K, 


where 


^C 


^c  = 


u. 


Y  = 


3h 

u 

K 

K2    = 

U 

Vh    = 

3C 

-F 

0 

G 

3h 

J 

t        . 

^3  n 

For  the  rest  of  this  chapter  we  will  be  dealing  only 
with  this  dimenslonless  moving  coordinate  system  unless 
otherwise  mentioned.   Thus,  from  now  on  we  shall  omit  writing 
the  circumflex  for  dimenslonless  variables  and  we  shall  use 
x,y,t  for  the  moving  coordinate  system  instead  of  5,  n  ,t. 

b .   Difference  Equations 

We  consider  the  rectangle  R  with  sides  of  length 
L-,  ,  Lp.   We  choose  a  rectangular  mesh  in  such  a  way  that 

the  coordinates  of  the  grid  are  defined  by 

IL 
X.  =  ^  i  =  0,l,...,I 

1    J  As  d     5  >    J 

where  the  boundary  lines  correspond  to  i  =  0,1;  j  =  OjJ. 

Li       L2 

For  the  unrefined  mesh  we  choose  I  =  -7 —  ,    J   =   t-     so  that 

As  '      As 

Xj  Y  are  integers. 

We  define  D.  as  the  connected  set  of  net  points 
in  the  interior  of  D  (i.e.  we  exclude  the  points  at  the 
northern  boundary  and  points  on  the  front).   By  a  regular 
point  we  mean  a  point  in  D.  whose  eight  nearest  neighbors 
are  all  in  D . ,  all  other  points  in  D.  are  called  irregular 
points.   At  regular  points  we  consider  two  different  second 
order  schemes.   The  first  is  a  one-step  scheme  that  is  a 
generalization  of  the  Lax-Wendroff  scheme  [12]  to  nonlinear 
equations.   The  second  method  is  a  two-step  scheme  due  to 
Burstein  [3]  . 
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In  order  to  formulate  the  one-step  method  we  write 
the  differential  equation  in  vector  form.   Let  W  be  the 
vector  of  dependent  variables  and  A  =  (a..),  B  =  (b..), 
and  C  =  (c^^)  be  matrices  that  could  depend  on  W. 


so 
where 


(c..  .) 


W,  =  AM^  +  BW   +  C 

W^^.=  A.  W   +  AM  .  +  B.  W   +  BW  .  +  C. 
tt    t  X     xt    t  y     yt    t 

Wxt=  ^W,  +  AW^^  .  B^Wy  -H  BWyy  -H  Cy 

\    =  ^^ij,t^'  h    =  ^^j,t^'   ^t  =  ^-lj,t 

Thus,  W  .  is  given  in  terms  of  space  derivatives  only. 
We  then  use  these  time  derivatives  to  calculate  W  at  the 
new  time  step  by  using  Taylor  series 

W(t+At)  =  W(t)  +  (At)W^(t)  +  (At)^W^^(t)  +  0((At)^) 

For  the  finite  difference  scheme  we  ignore  all  terms  of 
order  (At)   and  replace  all  space  derivatives  by  centered 
differences. 


Thus,  let 
^^f  =  f(x  .  _  _..,  ^,  .y. 


T^f  =  f(x  +  s  Ax,  y)      Tff  =  f(x,y  +  s  Ay) 


Then 


ax    2Ax  ^  X    y.  ''      9y     2Ay  ^  y    y  ^ 

2  2 

^    '     ■"■    (T  -2I+T"-'-)   -^  -y     — ^  (T  -2I+T"-'-) 


3x2     (^^j2   X     X    gy2     (^yj2   y     y 

2 
3^37  ^   r(Ax)(Ay)  ^"^x"  "^x  ^  ^"^y"  ^y  ^ 
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'   Let  X   =   max  (7—  ,  7—).   If  A,  B  and  C  were  constant 

Ax  *  Ay         ' 

and  if  we  were  to  consider  the  pure  initial  value  problem 

then  this  scheme  would  be  stable  if  X  <  ,  where  a  is 

~  a  /F 
the  largest  eigenvalue  of  either  A  or  B  (in  terms  of 

absolute  value)  [13]-   It  turns  out  empirically  (see  for 

example  reference  [2]) that  for  the  fluid  dynamic  equations 

this  condition  is  too  stringent  and  can  be  exceeded  in 

practice  without  causing  instability.   In  particular  it  was 

found  that  the  value  of  A  could  be  doubled  without  causing 

instability.   However,  in  using  this  criterion  it  must  be 

remembered  that  near  the  front  the  distance  between  the  front 

points  and  the  mesh  points  can  be  considerably  less  than  Ax. 

It  thus  seems  reasonable  to  modify  the  difference  scheme  at 

points  close  to  the  front.   The  eigenvalues  of  A  and  B  are 

u+c,  u,  u-c;  v+c,  V,  v+c.   Thus  we  require  that 

At  <  TT — I — I — n^ — ; — jr-        Where  As  is  the  distance  between 

-  max( |u I , I v| )  +  /h 

the  points  used  in  the  calculation  and  K  is  a  constant  to 
be  determined  by  trial  and  error  and  is  approximately  equal 
to  -7p— .   Near  the  front  where  As  is  small  h  is  also  small, 
so  that  there  is  some  compensation  for  the  short  distances 
involved. 

For  the  two-step  method  at  the  regular  points  we 
convert  the  system  of  differential  equations  to  conservative 
form.   A  partial  differential  equation  is  said  to  be  in 
cor  ■    ttive  form  if  it  can  be  written  as 

w,  +  f^  +  g,,  =  r       where  f  =  f(w,x,y,t),  g  =  gCw,x,y,t) 
u     x     y 
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Wit 


Ished  by  letting 


w  = 


hu 
hv 


f  = 


2   1   2i 

huv 
hu 


g  = 


huv 

V,  2^  1  .  2 

hv  +  p  h 

hv 


r  = 


Pu 

-Fv+G 
0 


A  two-step  Lax-Wendroff  scheme  Is  a  second  order 
;.eme  in  which  temporary  values  are  generated  by  a  first 
order  scheme.   Then,  in  a  second  step  this  intermediate 
data  is  used  to  generate  the  values  of  the  variables  at  the 
next  time  level  and  these  values  are  accurate  up  to  terms 
of  order  ((At)  ).   A  general  form  for  a  second  order  scheme  is 

.     .:    w(t+nk)  =  PjjW(t) 

:    w(t+k)   =  Q  w(t)  +  R  w(t+nk)  . 

These  two  steps  were  originated  by  Richtmyer 
as  generalizations  of  the  modified  Euler  method  for  ordinary 
differential  equations.   However  his  scheme  separates 
those  points  where  i+j  are  even  from  those  where  i+J  are  odd 
and  it  has  been  found  (see  for  example  reference  7)  that 
the  Richtmyer  scheme  can  be  weakly  unstable;  i.e.  there  can 
be  a  divergence  between  the  2Ax  and  ^Ax  components  of  the 
solution  because  of  this  lack  of  coupling  between  neighboring 
j.^lwvj.   We  shall  therefore  use  a  two-step  scheme  proposed  by 
in  [3].   In  this  case  n  =  1  In  our  general  formula 
ur  scheme  simulates  a  predictor-corrector  scheme  but 
with  only  one  '^^'^r'^ction. 
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"''^1+1/2 'yj  +  l/2>  *+"'  =  }   '"l,J  *   "l+l,j  *   "l,J+l  ^  "i+l,J+l' 


^   'f("l.lj   '  -  f^.J 


1  +  1, J+1  1,J+1 


At 


2Ay  'S^^i+l,j+l 


)  -  g(w 


•  _Ln   •    )  +  g(w.   .  ,  -,    )   -  g(w. 
1+1, J  &'   1,J+1  ^    1 


)} 


At 


^^^"l+l,j  +  l)  +  ^^"i+l,j-l)  -^  ^^^i-l,j  +  l)  +  ^^"i-i,j_i)> 


ww(x^,y  ,  t+At)  =   w(x^,y  ,t) 


At 


{f( 


w 


i  +  l,j   )  -  ^^"i-l,j   )  -^  ^^"i  +  l/2,j  +  l/2) 


-  ^^"i-l/2,j  +  l/2)  +  ^^"'i  +  l/2,j-l/2^  -  ^^"'i-l/2,j-l/2 


)} 


At 
^1A7  ^^^'"i,j  +  l 


{g(w.  ._^-,   )  - 


g(w.^j_l   )  +  g(wi+i/2,j+l/2) 


(w 


1+1/2, j-1/2)  ^    s(w._^/2,j+l/2)  -  s("i-l/2,j-l/2 


)} 


+  -^  {r(w.  .)   +   r.  .(t  +  At)}  , 


where 


r. . (t  +  At) 


IT  ^^^"i  +  l/2,j+l/2)  +  ^^"i  +  l/2,j-l/2^ 


+  Hw._-L/2,j  +  l/2)  -^  ^(«i-l/2,j-l/2 


)} 
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The  itrix  corresponding  to  the  linearized 


H 


(C.  n  )  =  I  +  IxIa  sin  U^^^L^^-IL)   +  b  sin  n  (3il^2^J.)l 

X2    r  1/2  1/2]^/^ 

+  —  .v'ACd-cos  C)(l+cos  Ti  )]    +  B[(l-cos  n)(l+cos  C)]    > 

The  eigenvalues  of  this  matrix  can  be  computed  numerically 

At    7550 
:  it  is  found  that  the  stability  requirement  Is  ^  <  '  ^ 

where  o is  the  maximum  sound  speed,  in  our  case 

o  =  /u^+v^  +  /h  . 

This  method  has  the  advantage  that  it  doesn't  involve 
matrix  multiplication  and  so  is  much  faster  and  it  also  allows 
a  slightly  larger  value  for  At.   However,  it  was  found  that 
::.     r  the  machine  time  was  spent  on  calculations  having  to 
do  with  points  on  the  front  and  so  the  time  saved  in  advancing 
the  solution  at  regular  net  points  is  not  large  compared  to 
the  total  time  required  for  the  solution.   The  disadvantages 
of  *:'-/'."    -ethod  are  that  it  is  not  a  disslpative  scheme  and 
that  it  does  not  generalize  to  the  irregular  points  where  all 
eight  neighbors  are  not  in  the  domain  of  interest. 

We  now  consider  ■'-•—  of  achieving  second  order  accuracy 
at  the  irregular  points.   Following  a  procedure  due  to 
Rlchtmyer  [l'<]  we  d '      the  irregular  mesh  points  into  two 
c^ci^ocs.   A  type  "A"  point  is  an  irregular  point  i  v- ' '^ '^  '■i^--Ar 
to  *     undary  than  a  certain  distance  6  (taken  to  be  one 
quarter  of  the  grid  distance).   A  type  "B"  point  is  an 
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irregular  point  not  of  type  "A".   At  type  "A"  points  the 
differential  equations  are  not  used  since  Gary  [I'l]  has 
found  Instabilities  If  this  Is  not  done.   The  reason  for 
this  Is  that  If  the  net  point  Is  very  close  to  the  front 
then  a  small  error  In  the  values  of  the  variables  on  the 
front  or  at  the  mesh  points,  will  cause  large  oscillations 
In  the  approximating  polynomial  based  on  these  values  and 
so  there  will  be  large  Inaccuracies  in  the  evaluation  of 
the  derivatives.   According  to  the  physical  interpretation 
of  the  stability  criteria  for  hyperbolic  equations,  based 
on   the  theory  of  Courant ,  Frledrichs  and  Lewy  [5]  one  would 
expect  trouble  whenever  the  distance  used  in  the  computation 
is  appreciably  less  thatn  c  At.   Thus,  if  we  would  use  the 
difference  equations  for  points  very  near  to  the  front  we 
would  have  to  reduce  At  to  an  unreasonably  small  value. 

To  avoid  this  difficulty  we  evaluate  u,v  and  h  at 
type  "A"  points  by  interpolation.  Instead  of  using  the 
finite  difference  equations.   We  first  find  u,  v  and  h  at 
all  regular  points,  at  type  "B"  points  and  at  the  points 
on  the  front  (by  methods  to  be  described  later)  .   In  the 
usual  case  we  first  interpolate  using  front  points  1,  2,  3,  ^ 
(see  Figure  ^)  to  find  u  and  v  on  the  front  along  the  same 
x-coordlnate  line  as  A, .   We  use  these  values  together  with 
those  at  the  mesh  points  B-j^ ,  R,  (but  not  other  type  "A" 
points)  on  the  same  x-coordlnate  line  as  A-,  to  find  u,v,h 
at  A-,  with  second  order  accuracy.   When  the  slope  of  the 
front  is  very  large  we  reverse  the  situation  and  do  the 

-29- 


FlGURh:  4 


u  el  A   is  ootained  by  interpolation  using  the  values  at  z  ,t)  , 
end  R  .  The  value  at  z   is  gotten  by  using  the  values  at  the 
front  points  1,^,^,A, 
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interpolation  along  a  horizontal  line.   In  all  cases 
encountered  here  it  was  possible  to  interpolate  in  either 
a  horizontal  or  a  vertical  direction. 

At  type  "B"  points  we  use  the  one-step  Lax-Wendroff 
scheme  independent  of  which  scheme  was  used  at  the  regular 
points.   The  only  change  is  that  now  we  must  use  noncentered 
formulas  to  achieve  second  order  accuracy  in  evaluating  the 
space  derivatives.   This  occurs  because  the  mesh  points  used 
in  the  calculations  are  no  longer  equidistant.   Thus  to  find 
u  at  the  point  Bp  (see  Figure  5)  we  first  find  the  location 
of  the  point  Zp  on  the  front  along  the  same  vertical  line 
as  Bp .   Generally  this  was  done  by  fitting  a  cubic  curve  to 
the  front  points.   However,  when  the  slope  of  the  front 
changes  radically,  as  when  the  occlusion  process  begins,  it 
was  found  that  a  cubic  polynomial  oscillated  too  much  and  a 
simple  linear  fit  gave  better  accuracy.   The  use  of  a  linear 
approximation  does  not  affect  the  order  of  accuracy  since 
it  is  used  at  only  a  few  points;  furthermore,  since  the  slope 
is  so  steep  the  front  is  nearly  a  straight  line.   In  all  these 
cases  the  same  number  of  points  were  used  on  either  side  of 
the  coordinate  line  to  prevent  distortions.   Thus,  at  most 
of  the  points  a  cubic  polynomial  was  used  even  though  a 
quadratic  would  have  sufficed  for  second  order  accuracy. 
Once  we  have  found  the  position  of  Zp  we  find  the  values  of 
u  and  V  at  Zp  by  interpolating  along  the  front.   We  then 
evaluate  u   at  B^  by  using  the  values  of  u  at  Rp,  Bp,  and  Zp. 
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u    is    obtained   at    B      by  usinj^   an  uncenterod    difference    scharae 
involving   R    ,fa    ,z    ,    u   9t   2    is   gotten    by    interpolation   along 
the    front  using    points    1,?,5,^^. 
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A  similar  process  is  used  for  the  x  derivatives  and  for  all 
the  second  order  derivatives.   Thus,  if  k  is  the  distance 
between  Bp  and  Zp  we  have  the  following  formulas. 


u 


B^  -    (Ay)UAy)  ^(^2)  ^  E^  ^^^2)  "  V^^  ^^^2^ 


(2.3)   u 


yy 


B2=  (Ay)(k+Ay)  ^^^2^  "  ITA^  ^^^2^  +  k(k+Ay)  ^^^2^ 


u 


xy 


u(A. )  -  u(R^) 

=  C ± 2 u 

Bp   ^      2  Ax         X 


B 


)  /  Ay 


2 


c .   Boundary  Conditions 

We  next  consider  the  finite  difference  approximations 
at  all  the  boundaries.   Along  the  east  and  west  boundaries 
we  use  the  regular  difference  approximations  using  the 
periodicity  condition  to  obtain  the  values  of  the  dependent 
variables  at  the  neighboring  points.   Along  the  northern 
boundary  we  have  v  =  0  for  all  time.   To  find  the  values  of 
u  and  h  we  use  one  sided  difference  approximations  and  arrive 
at  formulas  similar  to  those  in  (2.3). 

It  remains  to  satisfy  the  boundary  conditions  along 
the  front  as  given  by  equations  (2.2f).  Along  the  front  h  ==  0 
by  definition  and  so  we  must  solve  the  differential  equations 
for  u  and  v.   To  solve  this  system  of  first  order  ordinary 
differential  equations  two  different  methods  were  tried. 
The  first  method  was  the  trapezoidal  rule,  which  is  an 
implicit  method. 
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(2.'4) 


^.   +  At)  =  Xj^(t)  +  (At)<V^> 


V^(t  +  At)  =  V^(t)  +  (At)[Y<V>  -  <Vh^>  +  K2] 


Here  the  symbol  <  >  denotes  the  averaging  operator 

<w>  =  ^  (w(t)  +  w(t  +  At)),  and  y    and  K-  are  the  same  as 

in  equation  (2.2f).   Following  K.I.S.,  we  solve  (2.^b) 


for  Vj^(t  +  At)  with  the  result 


(2, 


■)  V^(t  +  At)  =  a  V^(t)  -  3<Vh^>  +  6  K2 


where 


a  = 


1+(|  At)^ 


1-(V^)2 


-F  At 


F  At 


i-(V^)2 


B  = 


1+(|  At)2 


At 


F(At)' 
w 


F(At) 


2 
At 


(c) 
(d) 
(e) 
(f) 


Set  Vh^(t  +  At)  =  Vh^(t) 
Predict  V^(t  +  At)  using  (3.3c) 


.:olve  for  the  dependent  variables  with  the  use  of  this 
implicit  scheme  the  following  procedure  was  used: 
(a) 

Move  the  front  •:--■"-  (3.3a) 

new  values  at  the  type  "A"  points 
.  :te  Vh»(t  +  At)  by  a  method  to  be 
If  there  is  any  -^  — '^^--^nt  change  in  Vp(t+At)  from 

ite  then  the  prO' 
star  b) . 
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Another  second  order  method  for  solving  ordinary 
differential  equations  is  the  Adams-Bashforth  method. 

Xj^(t  +  At)  =  X^(t)  +  ^[3V^(t)  -  Vj^(t-At)] 
(2.5)   V(t  +At)  =  V^(t)  +  ^[Y(3Vj_(t)  -  Vj^(t-At) 

-  (3Vh^(t)  -  Vhj^(t-At))  +  K2]  . 

This  method  has  the  advantage  that  it  is  explicit  and 

that  it  is  unconditionally  stable  while  the  trapezoidal 

rule  requires  iterations  and  is  weakly  unstable.   However, 

the  Adams-Bashforth  method  requires  knowledge  of  the  values 

of  the  variables  at  time  t-At.   At  time  At  we  can  avoid  this 

difficulty  by  finding  the  solution  by  a  finite  Taylor  series 

instead  of  the  finite  difference  scheme;  however,  when  we 

redistribute  the  points  along  the  front  we  no  longer  know 

the  values  of  u,  v  and  gradient  h  at  the  new  points  for 

previous  times.   It  was  found  that  until  the  time  that  the 

front  points  were  redistributed  for  the  first  time  both 

methods  gave  similar  results  and  so  it  was  decided  to  use 

the  implicit  method  only. 

Both  methods  for  integrating  the  ordinary  differential 

equations  require  the  evaluation  of  the  gradient  of  h  at  the 

front.   To  accomplish  this  we  first  evaluate  the  normal 

derivative  of  h  at  the  front.   The  slope  of  this  normal  is 

dy„  dy^ 

-1/  :3 —  ,  where   ^ —  is  found  by  using  a  quadratic  fit  along 
dx   '  dx  J  a         -i 

the  front  and  then  differentiating  this  polynomial.   We  draw 
a  straight  line  from  a  front  point  with  the  slope  of  the 
normal  into  the  cold  air  mass  (see  Figure  6).   We  find  where 
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The  valuos  of  h  at  7.  ,y.^   are  calculated  usinp^  the  values  at 
the  circled  points.  These  values  are  then  used  to  find^<^jL 


tMV 
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this  line  crosses  the  first  two  horizontal  coordinate 

lines  that  it  meets,  i.e.  at  z-,,z  .   If  the  first  coordinate 

Ax 
line  is  less  than  6  =  —r—   away  from  the  front  then  we  skip 

that  coordinate  line  and  consider  the  next  two.   We  do  this 
check  for  similar  reasons  to  our  not  using  the  difference 
equations  at  type  "A"  points,  i.e.  in  order  to  avoid 
Instabilities  caused  by  oscillations  in  the  approximating 
polynomials.   We  then  find  the  values  of  h  at  z-,  and  z„  by 
quadratic  interpolation  along  a  y-coordinate  axis.   If  the 
slope  of  the  front  is  too  large  we  again  switch  the  roles 
of  the  X  and  y  axes  and  find  the  value  of  h  at  the  inter- 
section of  the  normal  and  the  first  two  x-coordlnate  lines. 
Since  h  =  0  along  the  front  we  know  the  value  of  h  at  three 
distinct  points  and  so  we  can  find  the  value  of  the  derivative 
of  h  along  the  normal  to  the  front  with  second  order  accuracy. 
Also,  since  h  is  identically  zero  along  the  front,  by  defini- 
tion, the  tangential  derivative  along  the  front  is  zero. 
We  thus  know  both  the  normal  and  tangential  derivatives  of  h 
and  so  we  can  find  all  the  directional  derivatives  of  h  at 
the  front  and  in  particular  the  gradient  of  h.   In  fact  we  have 


The  sign  in  these  formulas  is  determined  by  the  direction 
of  the  normal  into  the  cold  air. 
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-i .       trlbutlon  of  the  Front  Points 

It  was  found  that  as  the  occlusion  process  continued 
the  individual  particles  on  the  front  converged  toward  the 
center  of  the  bulge  and  so  the  spacing  between  the  front 
--'nts  becomes  small  in  comparison  with  the  grid  size. 
This  leads  to  an  Inaccurate  exchange  of  information  between 
the  front  and  the  mesh  points.   In  addition  the  uneven 
spacing  of  the  front  point  caused  oscillations  in  the 

■oximating  polynomial  to  the  front.   It  was  thus  found 
necessary  to  redistribute  the  front  points  from  time  to 
time.   In  order  to  redistribute  the  particles  along  the  front 

at  time  t  we  calculated  the  arclength  between  particles.  This 

a   2 
was  done  by  passing  a  quadratic  curve  y-  =  p  x   +  3x  +  y 

through  the  points  i-1,  i,  i+1.   Then 


=   ax  +    6    ,  (;w^)      +    1   =   ax      +   bx   +   c 


dx  —    ■    M    ,  v^^ 

and 


^i"^i-l 


A+{^)^   dx   =    -^  ((2ax   +6)    /a 


^-1 


(j^)''   dx   =   ^-  ((2ax   +e)    /ax^+bx+c 

r-2 —  r^ 

+   log    |2ax   +   6   +    /ax   +bx+c|     ) 

i^-1 


When  the  ox^j,..  .;'  the  front  became  very  steep  we  considered 
^  f  y  and  arrived  at  similar  formulas. 

■/.trihuted  the  points  alont^  the  front  in 

---.-;.  ^  ... . .  ;..,..;.  Uw  .  ,  .;,/v..  were  equally  spaced  in  terms 

of  arclength.   Wf        he  new  values  of  x,  y,  u  and  v 
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by  interpolation  considering  these  variables  as  functions 
of  arclength.   Thus  everything  is  determined  within  an 
arbitrary  constant  which  can  be  fixed  by  specifying  that 
s  =  0  at  X  =  0. 

e .   Data  and  Results  of  the  Calculations 

The  following  numerical  values  for  the  parameters 
in  the  problem  were  taken. 

As  =  250,000  feet  ^  50  miles 

At  =  1,800  seconds^   p-  hour 

A   =  ^  =  .0072  second/feet 

X   =  20  As  or  about  1,000  miles 

Y   =  20  As 

C^  =  9.5  As 

C2  =  2  As 

g   =  32.1521  feet/second 

f   =  10   /second 


u   =  10  feet/second 

P_ 
P 


—  =  1  -  5I  :}:  .982 


u 


50 


P' 


feet/second 


So 


F   =  f  At  =  .If 


G   =  PA(^ 
P 


u'  -  u)  =  .05184 
=  3.1104  X  10"^  h 


and  in  case  I  we  have 


dh 

9y 


t  =  0    g(l-  ^) 


(p^  a- 


-    u)    '\j 


150 
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_    . .  aslty  discontinuity 

"   "    "-   -"ature  discontinuity  of  about  5  centi- 
grade, -.tationary  discontinuity  Is 

of  observed  values  for  the  slope  of  frontal 
surfaces  in  the  middle  latitudes.   The  time  step  was  chosen 

t  =  1/3  on  the  basis  of  the  stability  criterion  discussed 
i:.     previous  section.   This  time  step  was  reduced  by 
three-fourths  after  nine  hours  because  of  the  increase  in 
the  magnitude  of  the  velocity  In  the  domain  D. 

We  first  describe  the  result  for  case  I  where  the 
Ir'*-'-!  velocities  are  zero  in  the  moving  coordinate  system. 
F"      7  shows  the  initial  position  of  the  front,  which 
separates  the  domain  of  the  cold  air  from  that  of  the  warm 
air.   The  mesh  sizes,  in  the  unrefined  mesh,  were  taken  as 
h      length  and  width  1.   For  convenience  in  programming 
two  extra  columns  of  grid  points  were  added  outside  the 
western  and  eastern  boundaries  in  order  to  satisfy  the 

■ odlclty  condition.   All  the  diagrams  in  this  section 
will  refer  to  the        coordinate  system  except  where 
expressly  noted  otherwise.   Figure  8a  shows  the  position  of 
the  fr  ■:.      This  result  agrees  graphically 

v.-:  ■  K.I.S.       following  figures  show  the 

of  the  .;     -'*-   'in--  ii.^ur  intervals  until  a  total  of  16  hours 
.  tviu      ,       Mre  front  progress        e 

:Lnate  system.   The 
front  :        ♦•r..,.^  fa..u.  ,  than  the  warm  front.   Thus,  the 
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bulge  of  the  warm  air  Into  the  cold  air  narrows.   We  notice 
that  the  front  is  no  longer  a  single-valued  function  of  the 
X  coordinate  and  that  the  front  is  beginning  to  curl 
counterclockwise.   The  development  of  this  asymmetry  suggests 
the  occlusion  process  of  frontal  cyclones  which  agrees  with 
the  qualitative  analysis  given  by  Whitham  [19]  and  Stoker  [l6]. 
We  also  observe  that  after  approximately  l4  hours  a  second 
front  forms  to  the  west  of  the  original  front.   At  this  new 
front  the  warm  air  does  not  bulge  into  the  cold  air  as  far 
as  in  the  original  front.   Between  these  two  fronts  the  depth 
of  the  cold  air  is  quite  small.   Thus,  a  short  distance  above 
the  earth's  surface  these  two  fronts  merge  together. 

In  order  to  show  the  movement  of  the  front  in  detail 
the  trajectories  of  individual  points  on  the  front  during 
the  first  eight  hours  is  shown  in  Figure  9.   By  later  times 
the  points  on  the  front  have  been  redistributed  several  times 
and  so  it  would  be  difficult  to  follow  the  trajectory  of 
individual  points.   We  note  that  points  on  the  cold  front 
move  southeastward  and  those  on  the  warm  front  move  north- 
eastward on  the  average  whereas  both  fronts  propagate  eastward, 
The  movement  of  the  front  clearly  indicates  the  production 
of  a  cyclonic  circulation  about  the  circulation  center  where 
the  cold  and  warm  fronts  meet.   In  this  figure  we  also  show 
the  magnitude  of  the  velocity  components.   The  numerators 
are  the  x-component  of  the  velocity  while  the  denominators 
are  the  y-component  of  the  velocity.   Note  the  sharp  change 
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in  ■-  velocity  near  the  circulation 

■  '     ■■ .      To  Illustrate  the  circulation  pattern  even  more 
clearly  we  show  a  plot  of  the  velocity  field,  in  the  original 
coordinate  system,  in  Figure  10.   In  these  plots  the  direc- 
tion of  the  arrow  shows  the  direction  of  the  velocity  field 
while  the  length  of  the  arrow  is  proportional  to  the  magnitude 
of  the  velocity.   In  Figure  11  we  show  contour  lines  of  the 
height  of  the  cold  air  mass  at  5,000  foot  intervals.   In 
this  graph  Y  =  26  As  so  that  we  include  more  contour  lines 
and  so  that  this  corresponds  to  the  graph  of  K.I.S. 

The  trajectories  of  the  points  on  the  front  near 
the  circulation  center  shows  the  converging  motion  of  the 
cold  air  (see  Figure  9).   Consequently,  the  spacing  between 
these  points  on  the  front  becomes  small  compared  with  the 
grid  spacing.   This  uneven  spacing  Introduces  large  errors 
in  all  the  interpolation  formulas.   Therefore,  it  was 
necessary  to  redistribute  the  points  on  the  front  to  equalize 
the  spacing.   However,  it  was  found  that  if  the  redistribution 
was  done  too  often  the  results  were  smoothed  out  too  much 
and  the  deformation  of  the  front  no  longer  resembled  the 
occlusion  process.   Thus,  the  points  on  the  front  were 
redistributed  after  6  and  8  hours  and  then  every  hour 
afterwards.   Because  of  the  high  curvature  near  the  center 
of  circulation  it  was  necessary  to  approximate  the  front 
there  by  a  linear  function  rather  than  a  quadratic  or  cubic 
function.   The  higher  order  polynomials  had  excessive 
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oscillations  because  of  this  large  curvature.   For  the 
same  reason  the  velocity  components  near  the  circulation 
center  were  calculated  by  linear  interpolation  rather  than 
using  higher  order  polynomials. 

This  program  was  then  repeated  for  the  southern 
hemisphere.   Here  f  =  oa  sin  (f)  is  now  negative  since  (f)  is 
negative.   However,  as  initial  conditions  we  take  the  x 
component  of  the  cold  and  warm  air  velocities  as  negative. 
Thus,  initially  both  layers  are  moving  westward  instead  of 
eastward.   As  before  the  y  component  of  the  velocity  is 
zero  in  both  layers.   When  we  introduce  the  moving  coordinate 
system  it  is  now  moving  westward  instead  of  eastward.   In 
this  moving  coordinate  system  we  have 

F   =  f   At  =  -  F.,  ,    G   =  F   ^  (^  u'  -  u  )  =  G,,  . 
s     s         N'     s     sAsps     s      N 

As  expected  the  front  follows  a  similar  paytern  except 

that  the  occluded  front  moves  to  the  west  as  seen  in  Figure  12. 

As  a  check  on  the  program  and  also  to  improve  the 
accuracy  the  program  was  repeated  with  a  finer  mesh.   This 
can  be  done  in  two  ways.   We  can  keep  the  same  dimensionless 
moving  coordinate  system  described  in  the  beginning  of  this 
chapter  and  choose  the  length  and  width  of  the  mesh  sizes 
as  p-  thoir  original  value.   Then  according  to  the  stability 
criterion  we  must  halve  the  time  step  so  that  At  =  1/6. 
An  alternative  method  is  to  change  the  coordinate  system 
so  that  As  =  125,000  feet.  At  =900  seconds   and  X  =  .0072sec./ft 
and  keep   Ax  =  Ay  =  i.   Both  methods  were  tried  and  gave 
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t:i^  ^^..^    :...      -.  '^w^.y  obviously  should.  After  8  hours 

•3  results  t      -ally  identical  with  those  obtained 
^he  coarser  mesh  (see  Figure  13).   However,  as  the 
curvature  near  the  circulation  center  increases  the 

•Itional  front  points  give  greater  accuracy.   Similarly, 
in  the  cold  air  near  the  circulation  center  there  are  large 
gradients  in  the  height  which  are  measured  with  greater 
accuracy  by  the  finer  mesh. 

To  make  the  model  more  realistic  a  third  case  was 
i..--^^uced.   In  case  III  the  initial  conditions  were  changed 

■;hat  the  front  is  sinusoidal  only  in  the  middle  of  the 
region  of  numerical  integration.   Near  the  east  and  west 
boundaries  the  front  is  now  a  straight  line  parallel  to  the 
northern  boundary.   By  making  these  straight  portions  long 
enough  we  can  eliminate  the  effect  of  the  eastern  and 
western  boundaries   on  the  occlusion  process  for  the  times 
considered.   Thus  we  are  simulating  an  infinite  domain  in 
the  x-directlon  and  so  we  can  ignore  the  periodicity  condition 
which  was  artificially  introduced  for  mathematical  convenience, 
We  thus  have 

Case  III 

'C2  -  Cj_  0  <  X  <  K 

y^(x)  =<C^  -  C^  cos  (^^—^   x)      K  <  X  <  K+X 

C^  -  C^  K+X  <  X  <  2K+X 

rest  is  th'       as  in  case  I  in  the  enlarged  domain 
_  ■  _  K+X.   In  the  graphs  shown  here  we  used  K  =  X  =  20. 
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As  seen  in  Figure  Ih   the  occlusion  process  is  unchanged 
by  the  new  boundary  conditions.   However,  previously  the 
front  had  been  forced  northward  near  the  eastern  and 
western  boundaries  by  the  periodicity  condition.   Now  the 
front  remains  closer  to  its  initial  value  and  even  moves 
somewhat  southward  to  the  west  of  the  front.   This  movement 
Increases  the  length  of  the  cold  front  where  we  have  the 
steep  slope. 

In  case  II  we  change  the  initial  velocities  and 
height  distribution  to  conform  to  the  condition  of  a  geostro- 
phic  wind.   The  position  of  the  front  is  again  initially 
sinusoidal  and  since  the  initial  wind  field  is  geostrophic 
we  have  -rr-  =  0  and  -tt-   =  0.   Thus,  the  initial  slope  of  the 
height  of  the  cold  air,  -5—  ,  is  constant  with  respect  to  y 
but  is  variable  with  respect  to  x.   The  boundary  conditions 
are  the  same  as  in  our  original  formulation,  case  I.   These 
initial  conditions  are  physically  more  reasonable  than 
those  of  initially  constant  velocities.  However,  this  prolem 
was  more  difficult  to  handle  numerically.   Figure  15  shows 
the  position  of  the  front  after  11  hours  and  18  hours  have 
elapsed.   As  in  the  previous  cases  the  entire  frontal  system 
progresses  eastward  relative  to  the  moving  coordinate  system. 
The  cold  front  moves  faster  than  the  warm  front  and  we  again 
observe  the  beginnings  of  the  occlusion  process.   As  was  done 
previously  this  program  was  repeated  using  the  finer  mesh 


-^5- 


witl.  ^^.ow^.itially  i..-^    ^  -  •-  results  as  seen  in  Figure  16. 

his  case  our  calculations  do  not  completely  agree  with 
Kasahara  and  the  occlusion  process  is  not  observed  until 
much  later  times  than  is  noted  by  Kasahara. 
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f  Itnire    1  1i 


3    H01)R3 


?5.- 


-T — f — f — I — h~t — i — h~i—\ — f — I — < — i — f — f — i — I — I 

1  ^  5  7  9  1  1         13        15         17         19 

Initial   helfrht   contour  pattern   of   the   cold   air   for 
case    I.    The    contour  lines  are   drawn   at   5»003   foot 
Intervals.    Y   Is   equal    to    26^3   so   as   to   correspond 
to   the   graphs    In   K.I. 3. 
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figure  lib 


8  H0UR3 


Height  contour  pattern  of  the  cold  air  for  case  I. 
The  contour  lines  are  drawn  at  5fOOO  foot  Intervals. 
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figure    11c 


12   HCJRS 


?5     - 


or  . 


13  5  7  9  11         13        15        17        19 

Height   contor  pattern   of   the   cold    air    for  case    I. 

The   contour  lines   are   irawn   at   5,000    foot    Intervals. 
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flenire    lid 


16   HCWRS 


Helg-ht   contour  pattern   of   the   cold   air   for  case   I. 
The   contour  lines  are  drawn  at  5,000  foot   Intervals, 
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CHAPTER  III.   TWO  DIMENSIONAL  -  TWO  LAYER  THEORY 

In  this  chapter  we  discuss  problem  II,  as  given 
by  equations  (1.5).   As  in  the  previous  chapter  the  cold 
air  lies  above  a  region  D  of  the  x-y  plane  which  is 
contained  within  a  rectangle  R.   The  domain  D  is  bounded 
on  three  sides  by  straight  lines  and  on  the  fourth  side  by 
the  front  curve  C(t).   The  warm  air  lies  above  the  entire 
rectangle  R.   Thus,  in  the  domain  D  we  have  both  the  warm 
and  cold  air  masses.   The  variables  associated  with  the 
cold  air  are  denoted  by  unprimed  letters  while  the  warm  air 
variables  are  denoted  by  primed  letters. 

As  in  the  previous  chapter  we  consider  a  moving 
coordinate  system  and  introduce  dimensionless  variables, 
t 


I 

At 

? 

= 

x-ut 
As 

u 

- 

X  (u-u 

) 

X 

= 

At 

As 

n 

= 

y 

As 

V 

— 

Xv 

h  =  x^gci-  ^)h 

u'=  X(u'-u)  V  =  Xv' 

h'=  X  g(l-  ^)h' 

u  is  a  constant  while  At  and  As   are  units  of  time  and  length, 
Let 

.,  -    1     _    P 


F  =  f  At  ,    G  =  -  FXu  , 


1-  P-!-      P-P 


-79- 


We  then  use  (x,y,t)  Instead  of  (^,  n  ,t)   and  drop  the 
circumflexes.   Our  system  then  becomes 

u^    +   uu      +  vu     +   h      +    (r+l)h„   =   Fv 
t  X  y  XX 

I 
V.    +   uv      +   vv     +   h      +    (r+l)h      =-Fu  +   G 

t  X  y        y  y 


(3.1) 


^t  +  h^V^^  *  ""K  "■  ^^ 


u.    +u'v  +VU     +rh 

t                X  y             X 

'            '     '  II               ' 

V.    +uv  +VV      +rh 

t               X  y             y 


=   0 
=   Fv' 

=    -   Fu'    +   G 


(h'-h).    +    (h'-h)(u'+v')    +   u'(h'-h)^   +   v'(h'-h)„ 
1/  Ay  A  y 


=    0 


or  in  vector  form 


(3.2) 

where 


«t  "■  ^"x  +  BWy   =   f 


w  = 


u 

V 

h 

u' 

V' 

h' 

. 

f  = 


Fv 
-Fu+G 
0 
Fv' 

-Fu'+G 
0 


=  F 


V 

0 

-u 

1 

0 

+    G 

0 

V' 

0 

-u' 

1 

0 

0 

•                            ■ 

t                   < 

A  = 


u 

0 

1 

0 

0 

r+1 

0 

u 

0 

0 

0 

0 

h 

0 

u 

0 

0 

0 

0 

0 

0 

u' 

0 

r 

0 

0 

0 

0 

u' 

V 

h 

0 

u-u ' 

h'-h 

0 

u' 
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B   = 


V 

0 

0 

0 

0 

0 

0 

V 

1 

0 

0 

r+1 

0 

h 

V 

0 

0 

0 

0 

0 

0 

V' 

0 

0 

0 

0 

0 

0 

V' 

r 

0 

h 

v-v' 

0 

h-h' 

V' 

In  order  for  this  system  to  be  hyperbolic  it  is 
necessary  that  both  A  and  B  have  real  eigenvalues.   Since 
both  matrices  A  and  B  have  similar  structures  it  is 
sufficient  to  analyze  A. 

0  =  det  (A  -  XI) 

=  (u-A)(u'-A)  {-rh(u'-A)^+  [(u-X)^-  h][(u'-X)^-  r(h'-h)]} 


Thus  J  two  of  the  eigenvalues  are  X=   u,  X=  u'. 
For  the  other  eigenvalues  we  must  solve  a  fourth  order 
polynomial  equation 


[(u-X)^-  (r+l)h][(u'-X)^-  r(h'-h)] 


r  h(h'-h)  . 


Outside  the  region  D  there  is  no  cold  air  and  so  h  =  0. 
Then  this  equation  simplifies  and  we  can  solve  explicitly 
for  the  eigenvalues. 


X  =  u  (double)  , 


u  '  +  /rh  ' 


In  the  general  case  this  equation  can  be  solved  numerically 
It  has  been  found  in  all  the  cases  treated  that  the 
eigenvalues  are  real  when  u,u',  h,h'  are  real. 
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If  we  compare  this  system  to  that  obtained  when 
the  single  layer  model  is  considered  we  notice  several 
difficulties.   First,  this  system  is  no  longer  symmetrizable 
and  so  the  results  of  Lax  and  Wendroff  do  not  necessarily 
apply  even  to  the  linearized  equations.   Similarly  because 
of  the  interaction  between  the  warn  and  cold  air  masses  this 
system  can  no  longer  be  converted  to  conservative  form. 
Thuw  we  are  not  able  to  use  the  various  two  step  methods 
employed  earlier  but  one  must  use  one  of  the  more  involved 
one  step  techniques.   In  the  region  where  h  =  0  we  calculated 


that  the  sound  speed  is  c'  =  /rh '  .   For  the  parameters  used 

in  the  previous  chapter  r  is  approxomately  50.   If  at  the 

southern  boundary  h  is  only  twice  the  maximum  value  of  h 

(i.e.  the  maximum  height  of  the  warm  layer  is  twice  the  maximum 

of  the  cold  layer)  then  c'  [^  10c  where  c  Is  the  sound  speed 

of  the  single  layer  model.   Since  the  size  of  the  time  step 

is  inversely  proportional  to  the  sound  speed  we  must  use  time 

steps  that  are  about  1/10  as  large  as  those  in  the  single 

layer  theory.   This  together  with  the  necessity  for  one  step 

method  shows  that  even  with  modern  computers  the  time 

required  to  follow  the  front  for  a  reasonable  length  of  time 

is  quite  large.   As  an  estimate,  to  follow  the  front  for 

eight  hours  of  physical  time  with  a  coarse  20x20  mesh  would 

require  about  half  an  hour  of  computing  time  on  the  CDC  6600. 

For  a  'lOxiio  mesh  the  time  required  would  Jump  to  about  ^  hours. 

This  is  compared  to  about  5  minutes  of  computing  time  required 

to  follow  the  front  in  the  single  layer  model  with  the  finer  mesh, 
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According  to  our  original  assumption  the  warm  layer 
does  not  appreciably  affect  the  dynamics  of  the  cold  air. 
Thus ,  the  high  sound  speeds  in  the  warm  air  are  not 
physically  relevant.   So,  the  small  time  steps  are  made 
necessary  for  mathematical  rather  than  physical  reasons. 
However,  a  more  difficult  problem  arises  in  following  the 
motion  of  the  front.   In  the  single  layer  theory  the  front 
reduces  to  a  curve  in  the  horizontal  plane  and  can  be 
treated  as  a  free  boundary,  but  in  the  two  layer  theory 
the  front  is  a  contact  discontinuity.   Thus,  initially  the 
tangential  velocities  differ  on  the  two  sides  of  the  front, 
also  both  h  and  (h'-h)  have  discontinuous  tangents  across 
the  front.   This  discontinuity  of  the  first  derivatives 
probably  propagates  as  a  Jump  discontinuity  moving  through 
the  air.   It  is  well  known  that  along  contact  discontinuities 
Rayleigh  and  Taylor  instabilities  can  appear.   If  the  under- 
lying differential  equations  are  unstable  then  the  difference 
approximation  must  be  unstable  [15]  and  hence  these  physical 
instabilities  are  accentuated  by  the  errors  inherent  in  a 
numerical  approximation. 

A  first  attempt  to  solve  this  problem  was  made  using 
as  initial  conditions  the  assumptions  that  were  used  in 
constructing  the  single  layer  theory.   Thus,  we  assumed  that 
relative  to  the  moving  coordinate  system  v  =  v'  =  0,  u  =  0 
and  u'  =  u'  =  u.   It  was  found  that  Instabilities  occurred 
immediately.   Subsequently  it  was  found  that  the  situation 
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improved  when  the  initial  conditions  were  chosen  so  as 
to  satisfy  the  Jump  conditions.   Rlchtmyer  [1^]  also  found 
that  for  stability  it  was  necessary  to  insure  that  the 
initial  conditions  satisfied  the  Hugoniot  relationships. 
In  developing  the  single  layer  theory  we  used  the 
kinematic  conditions 

h.  +  u'h  +  v'h   =0,    h.  +uh   +vh   =0 
t       X      y  t      X      y 

or  after  subtracting  the  first  equation  from  the  second, 

(u-u')h^  +  (v-v')h^  =  0  . 
A  y 

"  w  be  the  velocity  vector  (u,v)  in  the  cold  air  and  w' 

the  corresponding  velocity  vector  in  the  warm  layer.   Then, 

this  equation  states  that  w-w'  is  perpendicular  to  Vh. 

Since  Vh  is  normal  to  the  front  this  equation  is  a 

restriction  only  on  the  normal  component  of  w-w'.   Thus, 

we  have  (w-w')   »—  =  0 .   Since  tt—  /  0  we  must  have  (w-w')  =  0, 
n  an  an  n    ' 

i.e.  the  normal  component  of  the  velocity  is  continuous 
across  the  front.   This  equation  has  meaning  only  as  we 
approach  the  front  from  Inside  the  domain  D  since  w  is  not 
defined  outside  of  the  domain  D. 

Let  y  =  yp(x)  be  the  initial  location  of  the  front. 
We  then  choose  as  our  initial  conditions  in  the  moving 
coordinate  system 

u=0  u'=u'-u 

V  =  (u-u')  ^         V'  =  0 

where   yc  '^  ^1  ■*■  ^2  ^^"  ^"To  ^^  * 
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With  these  new  initial  conditions  we  were  able  to  continue 
the  solution  for  several  time  steps.   However  instabilities 
still  occurred  after  about  six  time  steps  and  before 
meaningful  results  could  be  obtained.   Research  is  continuing 
to  find  methods  of  eliminating  these  instabilities  . 
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FORTRAN  I 

Two  DIML^ 

BURSTtlW 

COMMO.M/Al 

COMMOi^J/A? 

C0MM0.>J/A4 

C0MMni\J/A5 

C0MMGi\|/Af) 

COMHON/Al 

COMMON/Al 

COMMON/A] 

COMMON/Al 

COMMnN/Al 

CO^'MON/AJ. 

COMMON/A'f' 

COMMON/ A? 

COMMON/AP 

COMMON/A? 

COMMON/A? 

COMMON/A,'' 

CJMM0N/A3 

C0Mh0N/A3 

C0MH0N/A4 

COMMON/A^ 

COMHON/A'5 

COMMOiN/A"^ 

COMMON/PR 

CALL    CONI 

IF'     (IHLOT 

JPLOT    =    1 

K    =    1 

A  K  T    =    n , 

KT    =    Q 
AKT     lo    PiV 

KT  13  PKh 
CALL  ir:IT 
IF"  (IPLUT 
=    KT  +  3 


V    PR 

SION 

;iETh 

/AL, 

/AX, 

/F,G 

/U(? 

/N(2 

U/tR 

1/lP 

2/K 

b/L( 

b/lC 

9/OH 

1/Nl 

2/nF 

3/AK 

4/KT 

5/IP 

6/TT 

C/IT 

6/UI 

Q/FF 

Q  /  F  t^ 

i/FH 

2/FU 

1/LR 

;j  I T 

.UT 
APS( 


QlifrAH    FKOhUA(  INPUT, OLTF'iT) 
TlvO    LbVe.L    I.aV     WtrNLRCFF     .iFT-nD 

un 

AS,  AK 
AY 

4.27,3),V(24,?7,3) 
4,ki7,3) 

,  III,  JF,  Jrl 

24, ?7) 

hT,IpPl,IHK^,lPn,IPl,IP?,JPO 

HPK 

,  N  J  ,  (J  U  ,  ig  I  ^  ,  ivl  1  J 

,  NF1  ,MF2,riF.S,r  F4 

T 

HINT,  IPLOT 

, 'IT 

h'^,  ITFRT 

o  M  I  r  i 

(?b,6) 

x(on  ) ,  F  RY(oii )  ,F  RS  (>5n  ) 
X(3r ) ,FHY(iO) 
(  3  0  )  ,  F  V  (  0  n  ) 


,     0)     CAl L    PLTINIT 
IPLHT) 


fc<^tNT     TIMF     IN    NlWUTFS 
bEin    NL'IIFR    UF     Tl.ib    STtPS 

.UT,  C)  CALL  ."lYPt  OT 


KT 
IF" 
DT 
AK 
AL 
AS 


(KT  ,r 


=  ChPh: 
=  CHRP 
=      CFIPL; 

=  CHPF 
AKT  =  AKT 
IF  (KT  .F 
(\<T  .F- 
(Kf  ,F 
(KT  ,F 
=  KT-I 


h,  ICHT) 

K*UT 

R*AK 

K*AL 

K*AS 


no  TO  >5 


IF 
IF 
IF 
KZ 


J. 
CHT 


HT 

ICHT) 

ICHT) 

IPt^l) 

IPR2) 

*9 


IPKlwT  =  IPC 

jPLui  =  JPn 

IPkINT  =  IPI 
I  pw  I  I'll  =  IP2 
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20 
30 


35 


40 
50 


60 


70 


KL  =  IFTHtNtKT  ,r,J,     IChT.K/I.KT) 
IPRNOH  =  MOU(KL,  IPRINT) 
IPLNOW  =  dnU(KL,JPLOn 
PRINT  900, AKT 


IF 

DO 

IP 

IM 

DO 

IF 

L 

JP 

JM 

IF 

K 

WE: 

IF 

IF 


(IPRNOW  ,ECJ 
99  1=3, Ml 
=  I-l 

=  1-1 

99  J=1,NJ 


T)  PKIM  901,  KT 


(L( I , J) 

IS  kGUAL 
s  J*l 

=  J-1 
(K  ,GT, 
IS  bCUAL 


•  Lb 

TO  : 


0) 

IN 


GU  TO  99 
THt  Cni.0 


AIR,  0  IM  THP  WARM  AIR 


1) 

TO 


r,G  TO  20 
■?.    ON  THb 


CHECKING  IF  THt  NbAWfcST 
,ofc,  NJ)  GO  TO  99 


ARb 

(J 

(L( IP, JP)  , 


SbCHNO    TI'-'E    AROUNH 

IvblGhOORS    A''e 


IN    THF    COLD    AIR 


LF.n.OR 


L  (  I  P  ,  J  ) 


LE.O.nR.LH,  J°)  .LP.O)     GO    TO    99 


3  d  ,  t'  0  ,  7  U 
0) 


.LP 

.LF. 


CALL  TIMNbXKI.J) 
TIMNEX  IS  THE  TWO 
GO    TO    99 

( J-NJ*1) 

(L(  I'-),  JP) 

(L( I» JP) 

(L( IP, JP) 

IPUlNTsO 

IPOI.>JTsl 

IPUli>JT  =  ? 

IP0INTS3 

I  F  U  I  .N  T  =  4 

IPJIiMT  =  5 


STbP  8URSTEIN  MbThOD 


IF 
IF 
IF 
IF 
IF 
IF 
IF 
IF 
IF 
IF 


U(  IH,  J,>5) 
U(  Ifi,  J,>J) 
U(  I  • JM*3) 
U(  I  .  JM,«5) 


•  Lb 
THEN 

Then 
ThPr^ 

THEN 

ThF,\ 

ThFh 
I  POINT  =  IFTHfcN(L(  I,  Jrl)  ,Lb, 
IF  {L(IP,J)  ,LF,  0)  IpnlNT  s 
IF  (L(  If,  J)  ,GT,  0)  ('U  TO  4U 
IF  (  IPOINT  ,GT,  1  )  Gn  TU  tib 
IPOINT  =  2*IPun.T*3 
GO  'lO  6U 

IF  (L(  I'l,  J)  ,LF.  U  ,'JRi  L(  IM 
I  POINT  =  IFTHbf'd.  (  I,  JM)  ,L^-, 
PRINT  41?,  I.J,  IPOIM 
GO  TO  60 
IF 
CALL 


GO  TO  80 
0)  (-0  TO  61 

0)  GU  TO  '65 
NEI  THbR  J(I,wM,3),L(IP,J,"),U{I'^,J,M) 
U(  I  .  JM,3)  IS  Ml^^SINO 
IS  MISSING 
IS  MISSING 

Ann  L(IP,J,3)  Ai^E  HISSING 
ANL  L(  IM,  J,3)  A-'E  MISSING 
0,  1,0) 
2*IP0INT*? 


ARb    MISSING 


M) 


.LE. 
6) 


G)     3  0    TO    bP 


(  IPOlr.T.GT 


TIMNbXi;(  I 
TIMNEX  IS  THE 
GO  TO  99 
CALL  ONbSTEP{  I  , 
GO  TO  99 
CALL  NGT»-BON(  I  ) 


n,nf<,L(  IP,  JM)  ,LE  , 

.J) 

Twn    STf-p    biiRSTFIN 


0  ,lR,L(  IM,  JM)  ,LF,n )  GO  TO  60 


METHOD 


J. IPOInT) 
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80 
81 
82 
85 

86 
99 


105 

UO 


120 
125 

126 


130 


NOTHROM     IS    l-OH    LAUrULATIMG    U ,  V ,  t-     aT     THF    MjJ^THtPiM    BOUNDARY 

GO    TO    9V 

IF    (LP    ,t-iJ,     0)    Pf-IM    -^n^,     I  ,  J.HI  ,  J)  ,L(  I".  JP5 

L( I# J)    =    10 

GO    TO    9V 

IF     (LP     ,tnj,     0)     PRINT     4(16,     I,  J,L(I,  J),L(I,  J^M 

L(  I, J)     =    11 

GO    TO    9  9 

IF"     (LP     ,IU,     0)     PPirxiT    407,     I  ,  J,L(  I  ,  J)  ,L(  If".  JP5 

L(I, J)     =    1? 

GO    TO    99 

IK    (L(  I,  JM)     ,GT,     (1)    CO    TO    ^J6 

IF     (LP     ,f-y,     0)     PPIMT    41  U,     I,  J,  IpniNT,L(  I  ,  J)  .l-(  IP.  J)'L(  IM,  J) 

L(I,J)    =    9 

GO    TO    99 

IF     (LP     ,  1-  Q  ,     0  )     Pf-'lf^l    mi,     I  ,  J 

L(  I, J)     =    13 

COMTIrJUt 

IF     (K     ,GT,     1)     lU;     70    lu5 

K    =    2 

CALL    PF.^I0D(2) 

GO    TO    9 

K    =    1 

ITER    =    0 

ITfch  =  ITtR+1 

FRONT  AND  FKONAK  CALCULATE  THb  POSIT  lU^  0."  THk:  FPON'F  AT  TIME  I 

CALL  FRU^'T 

CALL  FPUr  AM 

DO  130  1=3, Nil 

DO  130  J=1,NJ 

CALL  POSINC  I  ,  J,  IPijS) 

IF  (  IPUb  ,GT,  n  )  (jO  Tij  120 

L  (  I  ,  J  )  =  0 

0     IS    FOR     THt    WARM    A  IP    AND    1      IS     FOR     THE    CO.P     A  IP 

U (I  , J , 3 )     =    0 , 

V  (  I  ,  J  ,  3  )     =    U  , 

H  (  I  ,  J  ,  3  )     =    U  , 

GO    TO    13n 

IF     (L(I,J)-1)     1.2'>, 1311,126 

L( I , J)     =    2 

WE    mPE    CALCULATIf'G    U.V,H    a"!     FCINTb     T^Q    inj^aP    Tu     Tt^E    FROMT 

CALL    TOuriEP(I,J) 

IF     (M(  I  ,  J,3)     ,';T,     0  ,  )     1>U    TO    130 

IF     (ITER     ,LE,     1)     PRiriT    700,     H(I,J,3),I,J 

H(  I, J,3)  =  ,00001 

CONTINUE 

DO  140  J=1,NJ 

L  (  1  ,  J  )  =  L  (  N  I  .  J  ) 


-91- 


140 


L(2,  J 

L(NI? 

L(N13 

CALL 

CALL 

tront 


143 

150 
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160 


165 


ir  (  I 
IF  (  I 
ITfcrtT 
IF  (  I 
ITEr^T 
GO  TO 
IF  (  I 
PKINT 
DO  1ft 
DO  16 
IF  (L 
CONTl 
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L(M2 
L(NI3 
IF  (L 
1  (A»\T 
IF  (  I 
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IF  (D 
CALL 
CALL 
IF  (L 
CALL 
CHGFR 
IF  (  I 

DO  ■^.'? 

DO  't2 
U(  I,  J 
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IF  (A 
CALL 
FO«MA 
1  1?H 

406  FOF^ilA 
1    llh 

407  FORMA 
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220 
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FROf 
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TEHT 
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TeK 

=  0 
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0  Is 
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»J) 
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^     .fc 
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PRHCI 
=  (F 
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KELA 
FROW 
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CHCiF 
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5     I 
3     J 
il) 
,1) 
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T(jr 
L(  If" 
r(3r 

L(I. 
T(jn 

L(  IP 
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=  L(i 
=  L(4 

on(  J) 

LCULA 

.Lt, 
.LE, 
I^   W 

,GT, 


TES  H  Sub  X  AND  Y 
TO  BE  UibU  AT  Thfc 
1)  uO  TO  145 

0)  GO  TO  150 
E  AkF  10  ITERATE  AGAIN 
10)  no  TU  15? 


AT  T^t  FRJK'T 

NEXT  FRONT  CALCOLATin.N 


.Lt,  4)  GO  TU  154 
,  ITtR 
3,  Mi 

1,NJ 
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1,NJ 

L  ( iN  I  ,  J  ) 

L  (  N  U  ,  J  ) 

s  L(3,J) 

s  L(4,J) 

0,  0  ,AtiD, 

60U,  ,UP,  AKT,tQ,72n,  .Ofi,  AKT  .  F  J  ,  9ftO  .  )  )  CALL  hYPRNTti(i) 

^    ,E0,  n)  CALL  MYPHNTl(.l) 

RS(NF3)-FRS(3) )/NF 

N  ,GT,  ,75*l)ELS)  f.O  TO  220 

8Lb 

FN 

Q,  0)  CALL  nyPHNTK?) 

R 

CKb  IF  TRX  IS  LESS  TI-AN  ZERO  OR  ii.^tATER  THAN  Nj 

w  (EC,  H)  CALL  MYPLOT 

=  1,M3 

=  i,;,j 

=  U( I, J,3) 

=  V( I , J,3) 

=  H  (  I  ,  J  ,  5  ) 

LT,  TT)  L.n  Ti)  1 

IT(1) 

H  tRK')R  HFSSA(it  1 

»JP)  =  ,I?/10(1H*)) 

H  tRROlt  MFSSAt.b  2 

JP)     =     ,  I2/1')(1H*)  ) 

H    bPH')R    MESSAUt    3 

,JP)     =     .I?/1U(1H«)) 


Al  POP'T  (  ,  I  ^,1H,  ,  13,611)  L  =  .\i, 
AI  POr-iT  (,I"^,1H,,I3,<5^)  L  =  ,\i., 
AT    pnr'T     (  ,  13,1  H,  ,  M,bH)     L    =     ,\i., 
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(  .  I  3  ,  1,  H 
L  (  I  M ,  J  ) 

!T  (,  n,iH,,  13, 

'A^^•1   i^iR  duT   pniMT 


410  FOPhATCJnH    tF^hOR    MF-SSAfit    4  AT     POT'T 
1=    ,I2,5H    L    =    ,Ik,llH    LdP.J)     =     iI2,llH 

411  FQRMAKvSflH    bRHOK    MFrSSAGt    5  AT     POI 
1    73H)        POINTS     ( IP, J) , ( IM, J)     UtPfc     IN 
2IN    COLD    AIR) 

415    rORi-lAT(iinx,llH    AT    POInT     (  ,  I  «i .  1  H  ,  ,  I  ?  , 

1    o4h)       L(IP#JP)     I'AS    f-'lSSlNf-       IPLIM    =    ,12) 

700  FORMAr(iX,*H  IS  STILL  LfeSS  THAN  ZfcPC  AND  =TUAL 
1  llh  AT  POINT  (,  I2,1H,  ,  I^^lfD/bUdH*  )  ) 

900  F0HMAT(27H  TML  MH-lPeH  UF"  MliMUTF.b  IS   ,F7,?) 

901  F0Pr,AT(^9H  THb  hI'MHfel'  0^'    TlnP  SIFPS  IS  ,lj) 
9in  FUPMAT(lX,*ITbi^  TS  (jK^AThH  THAN  3  AKH  FQUaL  TO 

END 


,I3,i4H)     HJInT 
=  ,I2/11(1H*)) 


(  I  ,  J  M )  WAS 


TO 


Fl?,7, 


I?) 


SUBRO 
COMI-iO 
CuMHO 
nOMMD 

c  u  M  r-i  n 

CUM  M  0 

CuMHO 

CUMHO 

COMhO 

COMHO 

COMMO 

CUMhO 

COMHO 

CUMHO 

COMMO 

CGMhO 

COMMO 

RtAU 

RfcAD 

data 

DATA 
DATA 
UaTa 
DATA 

data 
Data 

DAT  A 
Data 
DATA 
DATA 
DATA 


UTlr; 

N/Al 
N/A3 
.J/A4 
N/Al 
N/Al 

'Nl/Al 

N/A? 
i4/a2 
N/A2 
.-J/AP 
.M/A2 
N  /  A  3 
N/A3 
M/A3 
N/PL 

,j/pn 

bOU  , 
601, 

F,../ 
UT/l 
AX,  A 
TLtf' 
KDIS 
UISM 
CHPF- 
ICHT 
IPKl 
IPU, 
bPS, 


E  CO 
/AL. 
/AX, 

/F,n 

Q/tP 
d/IC 
9/CH 
0/TL 
1/IMl 
2/NF 

b/lP 
6/TT 
1/tP 

!3/Kn 

6/UI 
/SIZ 

1/LP 
IPR 
TT 

J,NF 

.la, 

u./ 

Y/i, 
/20. 
/.^5 

IN/5 
K/,7 
/5i»/ 

,  IPR 
IPl, 

bPSl 


N  I  N  T  T 
AS,  AK 
AY 

S 

hT,IPRl,IPK2,IP0,IP1,IF?,JP0 

PER 

tN 

,  N  J  ,  M  1 1 ,  i^J  I  iJ ,  N  I  3 

,NF1 ,NF^,NF3,nK4 

UINT, I  PLOT 

,nT 
bl 

IS 

bMiri 

EX,^I7t:Y,TLf-NY 
IMT,  IPLUT.LP 
/21,iil,^5/ 

.n:?j  64/ 

,1,/ 
/ 
/ 
,/ 

5/ 

2/110, l^rf/ 
lP2/tt,'*,i/ 

/.onom,  .uooom/ 
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DATA 
AK  = 
AL  = 
AS  = 

JPO  ; 


<S00 
601 


AK 

TT 

DT 

Nl 

NJ 

NF 

IF 

IPRI 

IPUU 

TLtN 

IF  A 

CuPfc 

ICHT 

IPRl 

IPO, 

JPO 

EPS 

EPSl 

SIZE 

TLEN 

Nil 

r.'I2 

NI3 

MFl 

NF2 

NF3 

rjF4 

RfcTu 

FORM 

FOR^. 

END 


1 
A 
A 

S 
S 
S 

s 
s 
s 

p 

NT 

T 
I 
P 


I 

IP 
lb 
IS 

I 


IZPX, 

./3, 

K/AX 

K/AY 

IFTHE 

DbLTA 

THE  T 

THE 

THE 

THE 

THE 

IS 


SIZEY,TLfc.'iY/7,d74,i,/S,i3,/ 


T 
T 
T 
T 

NC 


IS  T 
IS  TH 
S  THE 
01  NT 

IS  TM 
S  THE 
IPP2 
1»  IP? 
hoh 
TME 
S  THE 

S  THE 
IS  TH 
NUl 
NIf2 
NI*3 
Nf  ♦! 
NF*2 
NF*  j 
NF*4 


N(  IPL 

T   A 
OTaI. 
iMfc  S 
UTAL 
OTAL 
UTAL 
NZtRO 
HE  NU 
t  NUM 

LENb 

IS  LE 

INS 

t  PER 

TIME 
ARE  T 

ARE 
OFTEN 
bPkOh 

THE 

LfcNG 

t  LEN 


HT   ,'i 

V     IS 
TJMfc 
TF.P     I 
NUHUfc 
NUMBf- 
r'UMHE 
WF    ^- 
fiBFR 
PER    0 
TH    0^ 
SS    TH 
TEAD 
rEMTA 

(KT) 
Mfc  TI 

HOW  U 
WF  H 
ALLU 

ns  AL 
O.D. 

TH  UK 

HTH  I) 


T,  0,R,10 
DELTA  X 
IN  MINUTE 
;>i  MINOTFS 
R  OF  POIN 
K  OF  POIN 
R  OF  FRDN 
L  I  H  I  iM  A  T  F 
OF  TIHfc  S 
F  TIME  ST 

THE  X  AX 
AN  HDIS  T 
OF  THh  DI 
bfc  WE  HFC 

WHEN  WF 
Mb  STEPS 
F7EN  Wb  F 
LOT  AFTEh 
wED  IN  TF 
LOWED  BFT 
b ,  TO  MOV 

THE  AXIS 
F  THE  Y  C 


00  ) 
AY  IS 
S 


f^ELTA  Y 


X  nit^ECTION 
Y  niPECTION 


TS  IN  THE 

TS  IN  THE 

T  POINTS 

MlST  Cr  THE  PKINTUUT 

TtPS  EFTWEEv  PRINTOUTS 

EFS  BETWEE*^!  PLOTS 

IS 

0  THE  RJI-'NCftRY  WF  USE  I  NT6  =?PnL  A  F  I  ON 

FFENFf.TIAL  EQUATION'S 

IJCF  AK  WHE**!  WE  VIOLATE  STABILITY 

HEDLCE  AK  RiCAUSF  OF  STa=«ILITY  NEbDS 

WFEN  WF  HEniN  PRINT  IhG  .'-lORE  OFTcN 

PINT  AT  TIMES  ICHT, IFRl, TPR2 

TIME  ICHT 
E  SLLLTION  DF  rKX(yx)=Aj 
WtFN  ITERATES  IN  SOLVING 
E  ThF  FROf^T 

TO  PE  PLOTTED   ( IN  INCHFS) 
0CP:ir,4TF  X  I\|  TERMS  OF  nblTA  Y 


RN 
AT 
AT 


(FiO.2) 


SUBROUTINE  PLTIMT 
CALL  HLUTS(250,2'^HI  , '> 
CALL  PLUT(0, ,-11 , ,-0) 
CALL  PLOTd,  ,1,  ,-J) 
PbTuPN 
END 


l^VlUP  TLPhFL   ) 
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SUBROUTlNb  IMIT 
COMMOinI/A:^/AX,  AY 
CQMMnN/A4/F,G 

CaMhON/A5/lj(?4,t7,>"^),  V(24,?7,^<) 
C0MM0iM/A6/H(24,27,3) 
CllMH0N/A15/L{21,?7) 
COMIiO:\l/A?0/TLbN 
COMhOi>J/A?i/NI,NJ,NIi.Nl^#NM 
CUMMnN/A22/NF,Mf  1 ,  MF  ^  ,  MFv5 ,  M"  4 
COMMO,M/A40/FF  (?&,6) 
COMHON/A'50/FRA(3n)  ,  h  k y  (  io  ) ,  I- Kb  (  oO  ) 
CUf''HnivJ/A'51/FHX(3P),FHY(3n) 
CUMMOivl/Ab;:e/KJ(30  )  ,F  V(>5U  ) 
PIT    =    ?,*3,141^>V?65>i/TLtM 
DU    1     I=3,ivlF2 
AI     =     ( I-3)*AX 
FKX( I )     =     ,P*AI 

FKY(I)     =    9,:3*TL  brv?0.-.i*Tl..lr.M*f^C'^(PIT*FRy(  I  )  ) 
FHX(I)     =    -,l*TLfc?'«PIT*G*bir-:(PIT*rf<X(  I  )  ) 
FmY(  I  )     =     u 
FU(  I  )     =     0, 
F  V  (  I  )    =    (1  . 
CALL    PF:K2(1) 
PKS(2)     =    0, 
nu    2     1=3, NF3 

CALL    FROMDIS(  I , u , , 0 , , U 1 h  D I S  ) 
FkS(  I  )     =    FRS(  I-D+nit-  uIS 
CALL    Ptrt?(-1) 
CALL    FPuf'PO 

FRX     lb     rut:    X    COUPDINATl-    Ul-      [Hf     1-TH    POIMT 
FhX     li    il    bUti    X    AT    Thf     FKCjnT    PCINTb 
FU     IS    U    AT    THb    FFfjNT     PUINT^' 
DU    6     I =3, Mil 
■  ^-^  )  *  A  X 

'^.5*TLf;:N/2n,-,l*TLtN*C0b(MT*Al) 
Ul-     THh    FRONT    AT    (.njFT)I 


AI     =     (  !■ 

FF(  1)     = 

Fl-     IS    THf     FUSITiriH 

ti    J  =  l,r]J 

=    (J-1)*AY 

,  J,l)    =    0, 

, J,l )    =    U, 

(  AJ     .Lfc,     PI-  (  I  )  ) 


■J/ 


■aF    LiNtS 


nu 

AJ 

11(1 

V(I 
IF 
H(I 
L(I 

IF 


L^n 


TU    7 


, J,l )     =    b*( Aj-FF( I) ) 
,J)     =    1 

L    =     1     THbN     THI-    PUll/r     13     IM    TFP    CCI.)    AM 
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10 


20 


s  0, 

0 

ThbN 


GO  TO  8 

H( I, J,l) 

L( I, J)  = 
IF  L  =  0 
CONTINUE 
CALL  HErRIOPd) 
DO  10  Jsl.KJ 
Ld,  J)  =  L(NI,  J) 
L(2,J)  s  L(NIl,J) 
L(NI2iJ)  «  L(3,J) 
L(NI3,J)  s  L(4,J) 

1=1, M3 

J=1,NJ 

M  =  ?,3 

M)   =  U(  1  , 

V( I , J,l) 

M)  s  H{ I, J,l) 
CALL  fiYPRNTl(^) 
RETURN 
END 


ThF  POI.NT  is  In  THE  WA^H  AH 


DO  <iO 
DO  20 
DO  20 
U( I, J, 
V( I, J,M) 
H( I, J, 


J.l) 


SUBRO 
TIMht 
COMhO 
COMMO 
COMMO 
COHHO 
COMHO 
(  I,  J, 
H  (I  ,  J 
i-,5*A 
2'»-H(  IP 

4*h(  I, 

U(  I  ,  J 
1*H(  1, 

J-H( I , 
4-rl(   I, 

0*H(  IP 
7.,b«A 
ds.H(  IP 
V-h(  1  , 
A*,2'j« 


UTir- 
XI  I 
N/Al 
i\J/A4 
N/A5 
N/A6 
N/Al 
2)  P 
,2)» 
L*(H 
,J,i 
a*(H 
JP#1 
,2)s 
JP»1 
L*(H 
JP,1 
J,l) 
H{  IP 
,  J,l 

.  J*l 
J,l) 
AK*F 


fc  TI 
S  TH 
/AL, 

/r,G 

/L'(2 
/h(? 
1/IP 

fcPKE 

(IP. 
)*U( 

(IP# 
)*V( 
(.25 
)*u{ 
(IP. 

)*U( 
•U(  I 

#JP. 

)*H( 

(IPi 
)*U( 
•  U(  I 
•(M( 


I-NEV1(  I  ,  J) 
E    FIRST    isTcP 
AS,AK 


OF  THE  TWC-STEP  PURbTEIN  HETHOn 


4,27,3) ,V(24,?7, 3) 

4,27,3) 

,  IN, JP, Jn 

SPNTS  ( I  *     1/?, J  ♦  1/?,^) 

(H(1P,JH,])*H(IP,J,1)*P(I,JP.1)*H(I,J,1)) 


JP,1 
IP,  J 
JP,1 
1,  JP 
•  (  h  ( 

I,  JP 
JP,1 
I,  Jf^ 
,  J,l 
1)«M 

IP,  J 
JP,1 
IP,. I 


•  U(  IP,  JP,1)-H(  I  ,  JF,  !)♦')(  I,  JP.l) 
1  )-H{  I. J,1)«U(  I, J,l  )  ) 

•  V(IP,JP,l)-HnF,J,l)*V(lP,J,l) 
1 )-H(  I , J,1)«V(  I  , J,1  )  ) 

P,  JP,1)*U(  IF,wP,l)*.H(  IP,  J,1  )*U(  IP,  J,1  ) 
1  )*H(  I  ,  J,1).(J(  I,  J,l)  ) 
*U(  IP, JP,1)*U(  IF, jP,l ) 

l)*U{I.jP,l)*t-<IP,J,l)*J(IP,J,1)*U(ID,J,1) 
•U(  I  ,  J,l) 

ip,jp,i)-h(i,wP,i)*H(i,jp,n 

1  )  - h (  I  , J , 1) • H {  I , J ,  1  )  )  ) 
•U(  IP, JP,1)*V(  IF, JP,1  ) 
l)*V(IP,J,l)-^h(I,«P,T)*j{I,JP,l)«\/(I.Jo,n 

•  V (  I  , J , 1)  ) 


IP,JP,i)*V(lP,JP,l)*h(IP,J,l)»vMP,J.l) 


8+H(I,JP.l)*V(I.JP,l)*M(I,J,l)*V(I,J,1)))/H(I,J.2) 

VtI»J»2)=(,25*(h( IP, JP,1)«V( IF,wP,l)*H(IP,J,n*V(IP.J,l) 
1+H( I, JP,1)*V( I, JP,1)*H( I , J,i)*V( I, J,l) ) 
2-,5*AL*(H(IP,JP,l)*U(IH,jp,l)*V(IF,wP,t) 

3-H(I,JP,l)*U(I.JFM)*V(l.JP,l)*h(IP,J,l)*J(IP,J,l)*V(IP,J,l) 
4-ri(I, J,1)*U( I, J,1  )*V(  I. J.l)  ) 
S-,5*Ao*(M(IP,jP,1)*V(IH,JP,i)*V(TP,jP,l) 

6-H(lP,J,l)*V(IP,J,l)*v(lP,J,l)*h(I,jP,l)*^(I,JP,1.  )*^(I.JO,n 
7-H(I,J,i)*V(I,J,l)*V(i,J,l) 

8+,5*(H(  IP, JP,1)*M(  IP, JP,l)-h(  IP, J,1)*H(  IP.  1,1) 
y  +  h(I,JP,l)*H(I,JP,l)-H(I,J,l)«H(I,J,1. ))) 

A+,?5*AK*((.*(H(lP,JP,l)+h(IP,J,l)+h(I,JP,l)+H(l,J,i)) 
B-F*(H(IP,JP,l)*U(IP,JH,i)*H(IF,^,l)*ii(lP,j,l) 

CfH(  I,  JP,l)*iJ(  I  ,  JP,l)+h(  I,  J,l)»lJ(  I  ,  J,l)  )  )  WH(  I,,),?) 

RfcTURiN 


SUBRQUT 

TIMNPX2 

CUMMON/ 

CUMMOIM/ 

CUMMON/ 

COMMOi^l/ 

COMMON/ 

(  i,J,2) 

(  I  M  ,  J  ,  ? 

( I. JM,2 

(  IM,  Jil, 

H( I, J,3 

l-.2i?*AL 

2-^H(  I,  J, 

3+H( I , JM 

4-.25*AS 

t>*H(  I,  J, 

b+HC IM, J 

U( I, J,3 

1-,25*AL 

2-H(  IM, J 

3=H( IM, J 

4=H( IM, J 

t?+,3«(H( 

b+H( I , J, 

7+h( I , JM 

d-  .2t5*AS 

9=h( I , JM 


INE  T 
IS  T 
Al/AL 
A4/F, 
A5/IJ( 
A6/H{ 
All/I 

kfp 

)  KF-K 
)  KEP 
2)  Pb 
)  =  H 
*(H(  I 
2)*IJ( 
,?)*U 
*(H(  I 
2)*V( 
,?)*V 
)s(H( 
«(H(  I 
»1)*U 
,?)*U 
M ,  2  )  * 

IP,  J/ 

2)*H( 
*?)*H 
*  (  H  (  I 

,1  )*U 


IMMEX 
Ht  SF 
,AS,  A 
G 

24,27 
24,27 
P ,  H'  , 
REStN 
RhSEr^ 
RtSEM 
PRL-SF- 
(l.J, 
P,  J,l 
I,  J,^ 
(  I,  Ji^ 
,  JP,1 
I,  J, 2 
(  I  M  ,  J 
I,  J,1 
P,  J.l 
(  IM,  J 
(  I  M  ,  J 
II  (  IM. 
1)*H( 
I,  J.P 

(  I,  JM 
,  JP,] 
(  I  .  Jl 


2(I*J) 
COM'    3TI-P 


OF     IFF     TlvO-STF3    HMR^TFIM    METHOD 


,3) 

,3) 

JP, 

TS 

TS 

TS 

NTS 

1) 
*U 
-h 
2) 

♦  V 
-H 
2) 

♦  U 

♦  U 
1) 
2) 

JM, 
P, 
-H 
2) 

♦  U 
1) 


,  V  (  2  4  ,  2  7  ,  3  ) 


Ji'l 


(  I 
I  I 
(  1 
(I 


1/2,  J 
1/2,  J 
1/2,  J 
1/2,  J 


1/2,2) 
1/2,2) 
1/2,2) 
1/^,2) 


(  IP. J,l)-F(  IM, J,1)*U(  r', J.l  ) 

(  IM,  J,2)*L(  I'l,  J,2) 

-H(  IM,  JM,2)«(j(  IM,  J'-S'?)  ) 

(  I. JP,1)-F( 1 . JM,1  )*V(  I  . JM,1) 

(  1. JM,2)*V(  I  , JM,2) 

-H(  in, JM,2)*V(  IM, JM,P) ) 

(  I  , J,l) 

(  IP, J,i)*L(  IP, J,1  ) 

*U(  IM,  J,1)*F(  I  ,  J,2)*M(  I,  J,'')*U(  I  ,  J.2) 

♦U(IM,J,2)+F(I,JM,2)*J(I,JM,2)*U(I,JM,2) 

^)*0(  IM,  J|v,2) 

J.l)-H(  IM,  J,l)*h(  V\    1.1) 

(  IM, J,2)*F(  IM, J,2) 

-il(  IM,  JM,2)*H(  IM,  JM,0)  )  ) 

(  I, JP,1)*V(  I, JP.1  ) 

*V(  I,  JM,l)fF  (  I,  j,2)*lJ(  f  ,  J»2)*V  (  I  .  J.2) 
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C* 

D* 

fc* 

V 

1- 
!• 
3- 
4- 
&» 

7- 

b- 

V* 
A* 
B* 
C* 

U* 
h- 

(i* 
S 
I 

R 
750  F 
1 
2 
E 


n  (  i  » 
H(  IM 
,5.i 
.25. 

H(  IM 

(  I  ,  J 

.25* 
H(  IM 

H(  IM 

H(  IM 

,25. 
H(  1, 
H(  I, 
H(  IM 
,5.( 
H(  I, 
H(  IM 
,5.A 

.25. 
F.(H 

.25* 
M  (  P' 
B  s 
F  (S 
bTuP 
QRhA 
6H 

17HS 
ND 


JM,2 

»  JM, 
K.f  * 
{H(  I 
,  J#? 
,  3)« 
AL.( 
>  J>1 
,  J#2 

,  JM, 
AS.( 
JM,1 
JM,? 

»  JM| 
h(  I  , 

J. 2) 

I  J,? 
K«(G 
(H(  I 

(  I  »  J 
(H(  I 

,J,? 

AL.( 

d  .r 


).u( 
2).U 
(h(  I 
«  J»2 
).V( 
(^"(  I 

M(  IP 

).U( 

2).u 


H{ 1, JP, 


).V( 
)*V( 

2).v 

JP,1 

*^-(  I 

)*H( 
♦  (M{ 

»  J  #  2 
.1).U 
.  J#2 
).u( 
SGKT(U( 
fc.  .5) 


I  J 
IM 
J# 
•V 
*^/ 
J# 
J# 
M, 
*^/ 
I^ 


»  JH, 
1 
( 
J 


»  J 

*  J 
IH 

*h 

J# 
M, 
.J 

♦  H 
(  I 
*U 
M, 


1 
1 
J 
J 

1 

M 


•V 

.J 

?) 

•  V 

•  u 
1) 
?) 


1) 

2) 

JM, 

,  J 
-H 
2) 
1) 
.J 


*^(i»jM,2)*^-(iM,j,?)«j(i'%j,p)«v(r^,j,?) 

;^).V(  IM,  JK.i:  )  ) 

( I, j»i) 

.«f  )*H(  I  ,  JK,^  )«V(  1  ,  JM,  ?) 

♦H(  IM,  JM,2  ).V(  IK,  JM,?)  )  )  )/w(  I  ,  J,.^) 

'(  I  *  J*l) 

J{  IP,  J.D.V  (  IP,  J,l) 

.V(IM,J,l)*t.(I,j,2).')(I,J,P).v'(I.J,2) 
-.V(I^,J,2)*h(I,wM,?)*j(I,J'^,?>«y/(l,J'-,?) 
,  i').V(  I*^,  JK,i)  ) 
((  I, JP,1).V (  I  ,wP.l) 

i*v(l,jM,i)*h(i,j,2)*V(I,J,?).v(I.J,2) 
|.V(1,JM,2)>H(IM,J,2).;(IM,J,2)*V(IM,J,?) 
,2)*V(  IM,  J^,<:) 
IH,1)-H(  I,  vM,l).h(  I  ,  J'M) 
i(  I,  JM,2).K  1,  jM,2) 

-H(  IM, JM,2).H(  IK, JM,?)  )  ) 


T(ll 

u  = 
OPKO 


rt  AT  PCj 
,F9,5,6 
UTiNt  T 


(I,JM,2)*H(lM,J,.?)fh(IM,Jl,2))) 

.  J,l) 

(  I,  J.2)*H(  I,  JK,2).U(  1  ,Jf'\,  I) 

J,2)*H(IM,JM,2).U{Ih',JM,?)))))/H(I,J,3) 

l,J,<)*u(I,J,^)*V(I,j,,M.v/(I,J,3))*SQrtT(Hi(i,j,3)) 
PRINT  Ton,     I  ,  v.bM.lK  I,  J,7)  ,  V(  I  ,  J.>i)  .H(  T  ,  J,  3) 

INT  ( ,  I2,1H, ,  I?,1GH)   ST49  =  •F9,5, 
M   V  =  ,t-9,5,eH   H  =  ,F9.?,3  3X, 
iMNtX) 


SuBROuTIMt  ONfcSTFP( I, J,  IHOINT) 
ONbSTbP  IS  FOR  POINTS  NfcAR  TMfc  h PCNT 

WHEKF:  IJN*?YMNt-TRIC  niFFEhFNCES  ARF  f'E'^UIRFD 
ANU  UShS  AT  HOST  Tnf-  POINTS  U  P  ,  jf' )  ,  (  I  ,  IP  )  .  M  M  ,  JH  )  ,  (  I  »  ,  J  )  ,  (  I  ,  J  ) 

RkAL  Kl.Kg 

COMMO  >I/Aj  /AL,  AS,  ^k 

C0MMC^/A3/AX, AY 

CUMh0ii/A4/F,G 

CUMM0t^/A5/L(?4,27,3)  ,  V(^'»,?7,3  ) 

CUMM0U/A6/h(24,27,3) 

COMMON/All/IP, IK, JP, Jh 

COMMO  J/A15/L(2'<,27) 

C0MM0i^/A35/KI)IS 

COMMO'<I/A40/KF(?5,6) 

M     =     1 
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ion 


?on 


Ai  =  ( i-:^)* 

AJ  =  (J-D* 
DAX  =  1,/AX 
DAY  =  1,/AY 
ntLX  =  AV 
IF  IPOINTs- 
IF  IPUIWT=G 

IK   iPOiin  =  i 

IF  IPUlNTs;? 

IF  IPUINT=3 

IF  IP0Ii'JT  =  4 

IF  IF0INT=5 

IF  (IPOIfT 

Vlt  MPt  Tf-YI 
IJVMIS    CALCU 

no  TO  (IPO, 
USX   IS  u  bu 

USY  IS  U  SU 
llj  IS  U  SUR 
Si=CTIUN    100 

I 


USX 

vsx 
tisx 
uxx 
vxx 

HXX 
IF 


=  (  U  (  I  P 

=  (  V  (  I  P 

=  (  H  (  I  F 

=  (U(  IP 

=  (V(  IP 

=  (H(IP 
(  IPO  I  NT) 


SbCTIorj    ?Q0 
CALL    nFKH-ST 
CALL    UVI1ISX 
UVMI?    CALCU 
Kl    =    FX 
K2    =    AX 
HbLX    =    FX 
IF      (Ki     ,LT, 


AX 
AY 


1    THE 

ThFN 

thfn 

ThFivi 

THFN 

THFN 

TnFlj 

,LF, 

NG    TU 

LATFS 

«    X 
d    Y 

T 

lb    F 

,  f-       tv 
;  J,M) 

.  J ,  ri ) 
»  J ,  r-1 ) 

,  J,M) 

»  J  ,  ^1  ) 

/  J  .  f'  ) 

4on, 

IS  F 
(AI,A 
(I.J, 
LATFS 


r  WE 
fjFI 
U(I 
U(  I 
U(  I 
U(  I 
U(I 

1)  b 

FliNi 
U    A 

no,  2 


iiXT 

HtN 
-U(I 
-V(  I 
-  h  (  I 
'2,* 

-^.* 

-2.* 

400, 

op   w 

J,l, 

1,KH 
U     A 


HAV 
ThtK 
,  JM, 

",  J, 
,  J^|, 

U     TU 

V  Th 
ND  V 
I.  U,s5 
UXX 
I'YY 
IS  u 
htN 
IHUI 
",  J, 
f',  J, 
",  J, 
"(  I  , 
V(  I  , 
H  (I  , 
buO 

'■ItN 

Kb) 
.UFX 

hU    V 


F     ALL 
U(  I,. 
3)     IS 


vJ) 

>i) 
i) 
3) 


IS 
IS 

A  ML 
ANC 


8    OF     THE    UFl\PE'>T    iJFir.HROPS 

M,^^)  .L(  IP,  J.i")  »0(  IM,  J,M)     ARE    MlSSINij 

MISSIKT, 

M 1  s  s  I N  r, 

Ml  SSI  NT, 
L  (  I  P  ,  J  ,  :^ ) 
L  (I  M  ,  J  ,  .M 


A  P  E    M  T  s  s  I  r  ij 


ino 

t     'itAh 
QM     TF 

no  .2no 

IS    u 

IS    u    s 

Si'b    X 

'-It  I  The 

NT    =    0 

M  )  )  /  (  2 
M)  )/(2 
M  )  )  /  (  k 

J,H)4.L 

J  ,  f  )  +  V 
J  ,  '^  )  +  F 


EST    FRDijT    POINT     TU    COLUUM    AI 
F    FHONT 

,knu),    iPorjT 

SLH     XX  UXY     IS     U     SU^     XY 

UE     YY 

T         IJYT    IS    J    Siib    YT         UTT    IS    U    SJb    F 

W     (IP,  J)     WCi^     (TM,J)     IS    MISS  IMG 

,1 

,*AX) 

,*AX) 

.•AX) 

(  I  H  ,  J  ,  ri )  )  /  (  A  X  *  A  X  ) 

(  IM, J,M)  )/( AX«AX) 

(  I  M  ,  J  ,  M  )  )  /  (  4  X  *  A  X  ) 


(IP,  J)  IS  Missiijn     I,b   IPniMT  =  2,4 


, VF  X,F 
ON  TF 


X) 

F  FFONT 


30  n 


D^l 

DKf? 

nK3 

USX 

VSX 

HSX 

UXX 

VXX 

HXX 

GU  TO  (4PQ, 

SECTION  3cn 

CALL  'MEKFST 


1,/(K 
1,/(K 
1,/(K 
K2*UK 
K2*UK 
(K3  -K 

?,*(r 

?,*(n 


RDIS*mX)  uO  to  999 
1*(K1  +  K2)  ) 
1*K?) 

2*(K1-^K2)  ) 
l*UFX*(Kl-ts2)*Ut^2*L 
l♦\/FX■t■(Kl-^2)*DK2*\/ 
2)*nK?*H( I, J,M)-K1* 
Kl.llf-  )(,[)iK?_*[H  I  ,  J,M) 
Kl«VFX-nK2«V(I,.l,M) 
DK2*h{  I  , J,h)*UK3«H( 
4nu,4no,&no,bno,4no 

IS  F  OK  UHL-N  (  If,  J) 
(AI,AJ,-1,KH) 


(  l,J,l')-^l*J^-6*U(   IN,  J,M) 
(  I, J,N)-K1*DK3*V( IM,J,M) 
DK3*H( IM, J, A) 
+  LK3*L(  IM,  J.'l)  ) 
+  LK2*V(  r-S  J,fl)  ) 
Ih',  J,f^)  ) 
,  b  0  U  )  ,  I  PC  n  T 
iS  MISSINH     I ,e   IPO  I 


T  =  3.3 
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Call  uvMibx(  I ,  j,'i,i^^,uFx,\jf  x,f)i) 

Kl  s  AX 

K2  =  FX 

DtLX  s  FX 

IF  (K2  ,l.T,  WDIS-AX)  GO  TO  999 

DKl  =  1 ,/(Kl«(Kl*K?)  ) 

DK2  =  1,/(K1«K?) 

DK3  =  1,/(K2«(K1*K2) ) 

USX  =  K2*UK1*0( IP, J|M)*( 

VSX  =  K2*i}Kl*V(  IP,  J,N)*( 

HSX  =  K2*UK1«H(  IP, J,M)*( 

UXX  s  2,«^(DK1«U(  IP,  J,M)- 

VXX  =  2,<»(rKl*V(  IP,  J,M)- 

HXX  s  2,*  (nKl.«H(  IP,  J,M)- 

IF  (IPOIM-4)  40n,40l',tJ0 


SECTION  400  IS  FOR  WH£N 
400  USY  =  (U(  I, JP,M)-U(  I, JM, 
VSY  =  (V{ I, JP,M)-V( I, JH, 
MSY  =  {H(  I, JP,M)-H(  I, JM, 
UYY  =  (U(  I, JP,M)-2,*"(  I  » 
VYY  s  (V(I,JP,M)-2,*V(I, 
HYY  =  (H(  I  , JP,")-2,*'U  I, 
GO  TO  600 

StCTION  •SOO  lo  FOR  ^^HtN 
500  CALL  .>lER[;STY(  I  ,  AJ,-1,KH) 
CALL  UVHISY(  I , J,Kb,UH  Y, V 
Kl  s  AY 
K2  =  AJ-FF( I ) 
IF  (K2  .LT,  RDIS*AY)  CiO 
OKI  =  1,/(K1*(K1*K?)  ) 
DK2  =  1,/(K1*K?) 
DK3  =  l,/(K2*(t'l*K2)  ) 
USY  s  K2*DK1*U(  I  ,  JP,'')*( 
VSY  =  K2»UKl*v( I , JP,")*( 
HSY  =  K2fcJKl*H(  I,  JP,I1)*( 
UYY  =  2,*(I:M«1i(  I,  JP.il)- 
VYY  =  2,*(DK1*V(  I  , JP.M)- 
HYY  =  ?,*(DM«H(  I  ,  Jp.M)- 
IF  (DtLX-K2)  6ni},600,5bO 
biJO  DELX  s  K2 

600  IF     (IPOUT)     700,601,601 

601  IF  (  IPOP  T  ,Gfe,  6)  uO  Tu 
UXY  r  (  ,r.«{U(  IP,.JP,N)-U( 
VXY  =  (  ,«;•(  V(  IP,.JP,M)-V( 
HXY     =     (  ,3«(M(  IP, JP,h)-H( 

GO    TO    fton 

602  UXY  =  (Obr-,5«{U( IM, JP,M 
VXY     =     (VSY-,5«(V(IM,Jp,,i 


K1-K2)*LK2«L{I,J,1)-K1*UK3*UFX 

^l-K2)*LK2*V(I,J,<i)-Kl*DK3«vrx 

Kl-K2)«LK2*l-(  I  ,  J,  1) 

UK2»li(  I,  J,M)*jK3*jrx) 

1M2*V(  I,  J,M)"»-JK3*;Fx) 

l)K?*H(  I,  J,M)  ) 

0 

{I,JM)     IS    MCT    MISSIN^'     I.F,      IP0INTsn,2,3 

M  )  )  /  (  2  ,  •  A  Y  ) 

N  )  )  /  (  2  ,  *  A  Y  > 

t^)  )/(2  ,*AY> 

J,f^)*L(  I,  JM,M)  )/(  4V«AY) 

J,^)*V(  I  . JM,M)  )/( 4Y»AY) 

J,")+l-(  1  ,  JM,:')  )/(  AY*AY) 


(  I,  JM)     IS    MI  SSI  IMG 


t-  Y) 


I  ,t        IPnpiT    =     1,4,5 


TO    999 


K1-K2)* 
K1-(S2)* 
K1-K2)* 
UK2«U(  I 
I3K2*V(  I 
IJK2*H(  I 


b02 

I  M  ,  J  P  ,  f-1 
IM, Jr,M 

h-i,  JP,M 

)  - 1.'  (  II', 
)-V{  IN, 


L"K2*l(I,J,-l1-'<l*DK3*jrY 
LK2*v  (  I ,  J,  •))-Ki*DK3«vrr' 
CK2*h(  I,  J,  -I) 
, J,M)*  3K3*JFY) 
, J,M)*UK3*;FY) 
, J,M) ) 


)  )/AX-J?X) /AY 
)  )/AX-VSX)  /(VY 
)  )/AX-HSX) /ay 

^M,M)  )«nAY) »DAX 
^M,M)  )  #nAY) •DAX 
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700 


BOO 


HXY 
GO     T 
SfcCT 

UXY 
VXY 
HXY 

uu   = 
vv   = 

flH    = 

ur  = 
VT   = 

HT    = 

I'XT 

UYT 

VXT 

VYT 

tUT 

MYT 

UTT 

VTT 

HTT 

'J(  I, 

V(  I, 

H(I, 

IF     ( 

AU    = 

AV     = 

UL    = 

SB    = 

IF    ( 

RtTU 

H(  I, 

F'RIh 

L(I, 

RbTU 

1    OH 
■<i    IHh 
960    FOR;-i 
IfcP    A 

fnd 


=     (HSY-i'^^CHdM,  Jp,M)-H(  I 

0    flOO 

ION'    7un     I 


85n 


999 


95  0 


=     .^5 

=     .^5^ 

=     .^? 

U(  i, 

V(  I, 

H(  I, 

-  IJ  J  * 

-UU* 

-HH* 

r     -JS 

=   -o<; 

=  -US 

=  -US 

=  -  H  5 

=  -HS 

=  -UT 

=  -UT 

=       -riT 

J, 3) 

U,.'^) 
U,3) 
H(  I  ,  J 

AHS( 
AHS( 

fHbr.^ 

mK*( 

sd   ,r, 

R.M 

U,3) 
T    96r 
U)     = 

RiM 

A  r  ( 1 1 

u  = 

?URHr 

AFdl 
NU    to 


\&- 

•  (U( 

•  (  V( 

•  (h{ 

J,M) 
J,M) 

vJ ,  r-i ) 
USX- 
VFX  = 
( I  i  S  X 
X*US 
X*US 
X*VS 
Y*VS 
X*  (U 
Y*(U 

•  USX 
*VbX 
*(US 
=  U( 
=     V  ( 

=  h( 
,3) 
u(  I  , 
v'C  \, 
IF  (A 
UL  +  S 
fc.     . 

=  0, 
,  1. 
IPO  I 

rl      AT 

*FV. 

UTIM 
H  AT 
UAL 


b     FOR    WHtN!     ALL     8 

H  H  F  ;nJ  1  P  (j  I  N  T  =  - 1 
IP, JP,M)-U( IM, JP, 
1P,UP,M)-V(  Ihi,  JP, 
IP,  JP,!-!)-ti(  li'l.  JP, 


NtAKFST    ''tljHsnK^    ARF     PRhSFNT 

M)-L(IP.JH,'<1)+M(1M,JP,M))/(AX*AY) 
M)-V(I"»JM,-I)+V(IM,JM,M))/(AX*AY) 
M)-h(in,JH,i)+H(rH,jM.M))/(Ay*Ar') 


V  v/  *  I  I  S 
VV  *  V  s 
♦  VSY) 
X-UII* 
Y-UH* 
X-UH* 
X-Uii* 
S  X  *  \'  s 
S  X  f  V  s 
-  (J  L  ♦  U 
-IIU*V 

A>V'~^Y 

I  ,  U  ,  rl 
I  .  U  ,  M 
I  .  J  .  M 
,LT. 
JiO)  ) 
J, 3)  ) 

U  ,^T 
URT(H 
5)    PR 


Y-Hj 
Y-HS 

-UU* 

uxx- 

UXY- 
VXX- 
VaY- 
Y)-H 
Y  )  -  H 
XT-V 
XT-V 
)  -  ^^^ 

)-^AK 
)+AK 

5  +  aK 
0.) 


X+K  * 
Y-K* 
HSX- 
V3X* 
VSY* 
VbX* 
VbY« 
H  *  (  U 
H*(U 
T  *  U  S 
T*VS 
*(i'X 
*UT* 
*VT  + 
*hT  + 
UU     T 


VV 

UU  +  G 

VV*F 

LbY- 

UisY- 

VbY- 

VbY- 

XX-t-V 

XYfV 

Y-VV 

Y-VV 

T+VY 

,b*A 

.5*A 

.b*A 

n    ns 


SY 

VV*UXY-HXX*"*V'^X 

VV*LYY-HVY+"*VSY 

VV^VXY-HXY-.-^u'^X 

Vv*VYY-HYY-"n)SY 

XY)-IISX•HSX-UU*^IXX-V?X*H'■^Y-VV*^XY 

Yy)-USY*MSX-llU*HXY-VSY*h-^Y-VV*HT'Y 

*LYT-hXT*F*(/T 

*VYT-HYT-F*JT 

T)-LT*H3V-UJ*MXT-VT*HSY-VV*HYT 

K* AK*LTT 

K*AK*VTT 

K*AK*HTT 

n 


,     A  V  ,  A  n  ,  A  V  ) 

( I,  j.ci) )  )/nbLx 

IfjT     VbU,      I,.,bH,l)(  T,  J,3),VM,  J,3),H(  r,  J,3),UFLX 


J  ,  H  (  I  ,  U  .  3  ) 

i\i  T  *  1 

PCTfJT      (;i;i,lH,,I 

b  ,  6  H       V    =    ,1-9,5,6 
t    UtiESTHH) 

POINT     (,Ii?,lH,,I 
T(5     .F1^.7/1'IC(1H* 


?,1UH)  STAS    =     .F9.5, 

H       H    =  ,F9,3,dH    ntLX    =     ,P9.5,1&A. 

2,47H)  v\     AM    LFSS     THAW    l<=Rn     IN    'J''cST 
)  ) 


SUhKOUTiNt    NOTHL^rJ(  I  ) 

NOTHf-UN     IS    UStn    FC-P    PUli^TS    UN    Tl-F    NCPTHEPM    d'^UMDARY 
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AT  T 
CQMh 
COHH 
CUMM 
COMM 
COMM 

COMM 
J  = 
DAX 
DAY 
M  = 
UU  = 
VV  = 
HH  = 
USX 
HSX 

vsx 

USY 
HSY 
VSY 

uxx 

HXX 

VXX 

UYY 

HYY 

VYY 

UXY 

HXY 

VXY 

UT  = 

VT  = 

HT  = 

UXT 

UYT 

VXT 

VYT 

HXT 

HYT 

UTT 

VTT 

HIT 

IJ(  I* 
V(I, 
H(  I, 
RbTU 
END 


CN/A 
ON/A 

GN/A 
ON/A 
OW/A 
0:UA 

o;j/A 
*)j 
=    1, 
=   1. 
1 

U(  I 
Mi  I 
H(  I 

=  .^ 

=  ,'-J 

=  0, 

=  ,"> 

=  .i> 

=  .!? 

=  (U 

=   (h 

=  0, 
=  (U 
=  (  H 
=     (  V 

=  .;? 
=   .5 

=  - 1 

-UU 

-UU 

-HH 

s     -0 

=     -U 

=  -u 
=   -u 

=  -rl 
=  -H 
s  -U 
=  -U 
=  -H 
J, 3) 
J, 3) 
J, 3) 


rUKUARY    UE    Oi^t    0N£    bICE-L     L  I F  FEPE'.' :ES 

l/AL.A?;,  AK 

?/AX, AY 

4/F,G 

^/U(?4,i:7,.-^)  ,  V(k:4,?7,3  ) 

(^/H(?4,27,3 

31/1P# If, JP, 

?l/M,f!j,Nl 


/AX 
/AY 


) ,  v(k:4,?7,3 

) 

1 .  N 1 2  *  M  3 


,  J ,  ^  ) 

,  J ,  n  ) 

*(U(  IP, J,M) 

♦(!-(  IP,  J,M) 

*(Un  ,  J-?,f' 
»(H(  I  , J-?,M 

•(V(  I,J•7,^< 
(  IP* J,M)-2, 
(  IP» J,M)-2. 


-ii(  1  M,  J,  h  )  ) 
-H{ IM, J,H) ) 

')-4,*u( I, J^ 

i)-4,*H(  I,  Jl* 

')-4,*V(  I,  JK 

*U(  I  ,  J,  ii)*L 

♦H(  I  ,  J, '•)+!- 


♦LAX 
*LAX 


^  )  ♦  3  ,  *  ij  (  1  ,  J  ,  11 )  )  ♦  Ll  A  Y 
,^)*3.*H(I,J,.0)*UAY 
,K  )  )*LAY 

(  IM,  J,M)  )«n4X*»iJ 
(  IH, J,K) )»TAy**2 

(  I  ,  J-<:i,i-')-?,*U(  I»  Jf'.M)*(.(  I  ,  j,,-|)  )*JAY*»P 

(I,J-2,IO-?,*h(I,JM,M)f(-(I,^,i))*0AY**? 

(I,J  =  2,f'')-?,*V(IiJ^»M))«PAY**<? 

♦(U(  IP,  Jt.'D-iU  IP,  J'1,M)-L  (  IM,  J,M*  J(  rS.).1,M)  )*U\X*DAY 

*(H(IP,J,M)•H(IP,JM,M)-t-(IM,J,M)^^(I'<,JM,M))*OAX*DAY 

'5*(V(IP,JM,H)-V(lM,JM,M))«nAy*nAY 

usx-hsx*f*vv 

♦  VSX6H'5Y-r*Un  +  U 

*(USXfVSY)-Ut'*HSX 

SX*USX-Ul'*liXX-VbX«UiiY-HXX  +  F*VSX 

F;X*USY-UI'*UXY-VbY«UbY-HXY.»-F*V3Y 

SX*VSX-Ul;«VXX-VbX»VbY-MXY-F*ll,SX 

SY*VSX-Ul'«VXY-VbY*VbY-HYY-F«MSY 

SX*(USXfVSY)-HH*{UXX*VXY)-l)SV»HSX-tlu*HXX-VSX»H'^Y 

SY*(l)3X  +  VSY)-HH*(uXY*VYy)-llSY»MSX.IlU*HyY-V?Y*H9Y 

T«lJbX-lJLi*UXT-VT»USY-HXT-F«VT 

T«VbX-UU*VXT-VT*VSY-UYT-F*IJT 

T«(USX*V?Y)-'-'M*(UXT*VYT)-LT«HSX-UJ*HXT-VT*hbY 

=    L/(  1  ,  J,M)+AK*uT*,b*AK*AK*LT 

=     U, 

=     n{  I  ,  J,ri  )*AA*hT*  .b*AK«  AK*hTT 
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SUHROUTIMb  UVHIbY(  I  ,  J,Kb,Uf-  Y,  VTY  ) 

UVMISY    UALCULmTcS    U    Ah!)    V    ON     THt     FRC'T     I'l    T.HF     Y    niRttCTIONJ 

DIMtNSIUN    DiST(4) 

CQM|'iOiVA3/AX,AY 

ruMMO.-l/A40/t'  F<?fj,6) 

COMHn,x(/A'5Cl/t-RX(jn),H-!Y(sJ'^).KRS(;jO) 

coMHO.>i/A?2/KJi:^o )  ,rv(on ) 

CUMMniJ/A10n/X(4)  ,F1(4) 

C0MM0.i/Aiai/F^(4) 

Al    =    (I-3)*MX 

nsTd)   =   0, 

DO    m     IK    =    ^,4 

IS    =    KU+IK-1 

DS    =    SQKK  (h  KX(  I-^-D-hRXC  IS)  )«(t-RX(  IS-l)-.-wx(  IS)  ) 
1    *(KRY(IS-l)=FKY(IS))«(PrJY(lS-l)-FRY(lSn) 
in    nsT(IK)     =    nbT(  IK-D-^IJo 

nSTX    =     JST(;^)*55(.'f>T(  (I- i<X(KLf^l)-Al  )«(F^X(Knf1.)-Al  ) 
1     ♦(FRY(KR  +  l)-rF(  I  )  )  »()■  KY(K<  +  l)-l-r(  I  )  )  ) 

no    ^0    KL=1,4 

Ki    =    KH+KL-i 

X(KL)     =     PSTCKL) 

FKKL)     =    FU(KS) 
20     F2(KL)     =    FV(KS) 

CALL     INTtKCUSTX,'','^,"!-  Y,  VKY,DtM) 


UFY  IS 
RfcTuR.J 
END 


U  UN  ThF  front  AT  X  COCRLIKATF  AI 


SUHROUT 

UVMISX 

D I  ML MS  I 

CUMMOM/ 

CuMMQiJ/ 

COMMON/ 

COMHOiv)/ 

COMMON/ 

COMMON/ 

COMMON/ 

COMMON/ 

Al  =  (  I 

AJ  =  (J 

CALL 

KA  = 

KK  = 


FD 
IF 
KB 


DbT(l) 


INb 

CALC 

UN  r 

A3/A 

AlO/ 

a5o/ 

Aion 

Ml  01, 

Kl/X 

f;?/m 

-3)* 

-1)* 
IS(  I 

TRtN 
+  KA 

=  0. 


UVH 

ULA 

ST( 

A  ,  A 

l:RS 

I-Ra 

MJ( 

/X( 

/F^ 

VPi 

F 

AX 

AY 

,J, 

(Mh 


IbXC  I  ,  J,  IDlKiKrf,Uh  X,VFY,FX) 
TbS    U     Aid)    V     01,^     THfc     FRCriT,     I  \l 
4) 
Y 

(^n  )  ,K^:Y(s5tl  )  »!•  RS(  JO  ) 
?  0  )  ,  F  V  (  Ji  (I  ) 
4  )  ,  F  H  4  ) 
(4) 


THE    X    UIKPCTinM 


IDTH,KH,FX) 
.l.b,     ^,0,1) 
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DO  10  IK«2,4 

IS  =  KP^IK-1 

DS  -    SQHT(  {^  RX(  IS-1)-FHX(  IS)  )«(FPX(  l'^-l)-"t'X(  I*^)  ) 
1  ♦(FRY(lS-l)-KRY(IS))«(rKY(lS-l)-fRY(IS))) 
10  DST(IK)  «  DbT(  1K-1)*I3S 

DSTX  s  UST(KA*l)*snKT(  (FH'X(KK)-XVF1)*(FRX(KK)-V\/P1) 
1  ♦(FRY(KK)-AJ)*(FRY(KK)-AJ) ) 

DO  20  KLsl,MF 

KSj  s  KB*KL-1 

X(KL)  =  CSTIKL) 

Fl(KL)  =  FU(KS) 
20  F2(KL)  =  FV(KS) 

CALL  lNTFrt(USTX,MF,2.UFX,VM,CiJf) 

RtTuP.J 

END 


SUBHOUTiNt:  FDIS(I,J,IUIR»KH,rX) 

FDIS  CALTULATtS  THE  DISTANCE  TO  ThF  FRONT  IN  TmE  X  JlK^CTION 

C0Mh0i>l/A3/AX,  AY 

COMMON/AlO/bPb 

C0^'M0i«J/A22/NF,  NF1,NF2,NF>J,MF4 

COMllOiM/A30/lTbR,  ITFNT 

COMNO,^^/A5  0/^  RX(3n),Fi'Y(in),FPS(iO) 

CUMKnN/AlOn/X(4) ,F1(4) 

COMHO,N/Fl/XVPi 

C0MM0N/F2/MF 

CUMMON/HPl/LP 


liO    TO    2 


1)     KP    =    fJFl 


AI     =     (I-3)*AX 
AJ     =     (J-1)*AY 
IF     (Kd     ,GT,     0) 
KB    =    1 
no    TO    4 
2     IF     (KB     ,RT,     NF 

4  ISPbC    3    0 
L    =    Kd*i 

CALL  JfcK(L,UERlV) 
IF  ( AdSCrtRlV)  ,nj 
LL  =  L*l 

CALL  JEK(LL»ntRlV) 
IF  (AdStDtRlV)  ,LE 

5  KB  s  Kb*l 
MF  =  2 

GO  TO  10 

6  MF  =  4 

Kb  IS  THF  FKONT  POINT 


1.7)  GO  TC  5 


1,7)  GO  TC  6 


WITH  WHICH  WE  ngniN  OUP  INTERPOLATION 
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10 


14 


60 


80 
86 


89 


140 


165 


200 


bQO 


DU 
X( 
Fl 
IC 
XV 
CA 
IF 
IF 
XV 
CA 
IC 
IF 
XV 

IF 

XV 

XV 

FV 

GO 

K 

IF 

K 

IF 

K 

KK 

KK 

IF 
IF 
GU 
IF 
KB 
IF 
KB 
IF 
KB 
GO 
IF 
IS 
IF 
GO 
FX 
Rt 
FO 


1  * 

600  FO 

1  I 

EN 


14  I 
INH) 
(  INU) 
OUMT 
Ml  = 
LL  IN 

IDIR 

iniR 

=  TH 
UL  IN 
OUNT 

<  ICO 
PI  = 

ARE 

(ABS 
Ml  = 

=  XV 
Ml  = 

TO  6 
=  KB 

(FKX 
=  Kfl 

(K  , 
=  KB* 

=  K- 

15  T 

(MF 
(ITt 
TO  2 
(KK 
T  =  K 
(ITE 
s  KB 
(KB 

=  1 

TO  ?. 
(KB 

PEC  = 

(  ISP 

TO  1 

=  A8 

TURN 

RilATC 

KB  C 

RMAT( 

DIR  = 

D 


1*MF 

RXCKR+IND-I) 

FKY(KB+IMU-1) 


AJ 


(XVM1,MF,1,FVM1,UUM,DCM) 
1  NE  GO  TO  THR  RIGHT 
-1  WE  no  Tf)  THE  LEFT 
F(IDIH  .GE,  0, , Al*AX, AI-AX) 
(XV,MF,l,FVX,DUM,DlJf  ) 
CUUNT+1 

,GT,  It?)  'iO  TO  80 
FVX*(XV-XVM1)/(FVX-FVM) 

NO  THE  HETHQU  OF  FALSE  POSITION  TO  SOLVE  FRX(XX)-AJ=0 
-XVPl)  ,LT,  EPS)  GO  TC  80 


i\in  = 

=  0 

AI 

TFK 


EMI 
TER 
'  I 
UMT 
XV- 
USI 
(XV 
XV 

PI 

FVX 
0 


(K)  ,GE,  XVPl)  GO  TO  89 

LF,  KB*3)  GO  TO  86 

4 

1 

HE  FRONT  POINT  SUCH  THAT  FRX(KK,3)  IS  LESS  THAN  XV 

AND  FRX  IS  AS  LARGE  AS  POSSIBLE 
,GT,  2)  GO  TO  140 

R  ,Lt:,  1  .AND,  LP  .EQ,  0)  PRINT  600i  I  #  J  i  I  D  I  Ri  K^»  X  V  ,  F  VX 
00 

,eu,  KBfl)  GO  TO  200 
K»l 

R  ,LE,  1  .AND.  LP  ,EQ,  0)  PRJNT  500i  1 , J# I D I R, KR, KBT, X V , F Vx 
T 
,GT,  0)  GO  TO  165 

00 

,GT,  NFl)  KP  =  NFl 

ISPEC+1 
EC  ,GT,  0)  MF  =  2 
0 
SF (Al-XVPl) 

IIH  AT  POINT  (, I2,1H,, I2,9H)  IDIR  s  ,12, 

HANGED  FROM  ♦»I2,4H  TO  ,I2,6H  XV  s  ,F11.7,7H  FVX  =  ,Fll,7) 
IIH  AT  POINT  (, I2.1H,, I2,39H)  WE  USED  LINFAR  INTERPOLATION 
,I2i6H  KB  =  ,I2,6H  XV  =  «F11,7#7H  FVX  =  ,F11,7) 
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SUBROUTIN 
SUBROUTIN 


10 


COH 
COH 
COM 
COM 
COM 
COM 
COM 
COM 
COM 
COM 
COM 
DIM 
FRX 
FHX 

DO 
IF 

OFU 
OFV 
OFR 
OFR 
OFM 
CFH 
3  FU( 

1  •. 

2  /( 

FV( 

1  ♦. 

2  /( 

TRF 

TRF 

FRX 

FRY 

IF 

CON 

CAL 

FRS 

DIS 

DO 

CAL 

IF 


MQN/Al 
MCN/A4 
M0N/A2 
M0N/A2 
M0N/A5 
M0N/A3 
hON/AB 
MOiN/AJ 
M0N/A5 
MON/A? 
MCN/Al 
ENSION 
IS  TH 
IS  H 
IS  U  A 
10  Is3 
(ITER 
(I)  » 
(I)  » 
X(l  ) 
Y(  I  ) 
X(  l) 
Y(  I  > 

I)  « 

5*AK*| 

1,*.25 

I)  =»  « 

25«F«A 

1,*.25 

X  «  FR 
Y  «  FR 
(I)  « 

(1)  * 

(AdS(F 
TINUE 
L  REH2 

(2)  = 
MIN  = 
100  I* 
L  FRON 
II.GT, 


£  FRO 

E  FRO 

A 

/AL,  A 

/F,G 

2/NF, 

4/KT 

0/lTE 

1/EPS 

6/DIS 

0/FRX 

1/FHX 

2/FU( 

00/X( 

CFU( 

E  X  C 

SUB  X 

T  THE 

,NF2 

,GT, 

FUU) 

FV(  I) 

FRX( 

FRY( 

FHX( 

FHY( 

(li-, 
FhX(  I 
«F*F* 

•  F*AK 
K*AK* 

•  F«F« 
X(  I  ) 
Y(  I) 
OFHX( 
OFHY{ 

Rx(  n 


NT 

NT  IS  TO  FIND 
ND  THE  VALUES 
S,  AK 


\F1,NF2,NF3,NF4 


TKE  NEW  POSITION 
OF  L  AND  V  Th:re 


OF  THE  FRONT 


R,  I 
1 

MIN 
(30 
(30 
30) 
4), 
30) 
OUR 
AT 
FP 


TERT 


),FHY(30),FRS(30) 

)  ,FHY(30) 

,FV(30) 

Fie*) 

, OFV(30),OFRX(30>,CrRY(3D),OFHX(30),OrHY{30) 

DINATE  OF  THE  l-TH  POINT 

THE  FRONT  PCINTS 
ONT  POINTS 


1)  GO  TO  3 


) 
) 
) 
) 

5*F*F*AK*AK)*0FU(I)*F*AK*0fV(!) 

♦  nFHX(I))-',25*F*AK«AK*(FHY(!)*0FHY(l).2,*G)) 

K*AK) 

0FU(l)*(l,»,25*F*F*AK*AK)*Drv(I) 

FHX(l)*0FHX(I))-,5«AK*(FHY(I)*0FHY(l)-2,*G)) 

K*AK) 


I) 
I) 
I) 
I) 

25 
)♦ 

AK 
*0 
(F 
AK 


1 )*,5*AK*(FU( I )«OFU(  I  )  ) 
I )*,5*AK*(FV( I )*OFV(  I  )  ) 
-TRFX),GE,&PS1,CR,ABS(FRY(I)bTRFY),GF,EPS1) 


ITERTxl 


(I) 

0. 
5, 

3,NF3 

DIS( I,0.,O,,DIFDIS) 

3  ,AND,  DIFIJIS.LT.DISMIN)  CISMI'yl  a  DIFUIS 
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100  FKSd)  =  FRSC  I-1)*DIKUIS 
CALL  Pf:R2(-l) 
RETURN 
END 


SUBROUT 
SUBROUT 
bETW 
COMMON/ 
COMMON/ 
SRP(XX) 
DSTINK 
1*AL0G(A 
ITIMES 
J  =  lAR 
IF  (I  , 
DX  =  PR 


rxo  = 

rxi  s 

rx2  = 

FYO  = 

FYl  = 

FY2  = 
GO  TO 
DX  = 

FXO  = 

FXl  = 

FX2  = 

FYO  = 

FYl  = 

FY2  = 
IF 


INE  ("R 
iNt  f  R 
EEN  TH 
A30/IT 
A50/FR 
s  SQR 
XX)  •• 
dS(SHP 
=  1 
SJI) 
LT,  0) 
X<J) 
RX(J-1 
KX(J) 
RX(  J*-l 
KY(  J-'l 
RY(J) 
RY( J*l 


ARCLENRTM 


GO  TO  1 


) 


) 


X 
=  FRX(J*1) 
s  X 

=  F 
s  F 
=  Y 

=  FKY(J+2) 
(FXO  ,NE,  FXl) 


•RX(  J  +  2) 
■RY(J*1) 


FXO  =  F 
FYO  s  F 
IF  (FXl 
FX2  =  F 

FY2  =  F 
0X01  3 
DX02  a 
DX12  = 

TERMO  = 
TERhl  = 
TERM2  = 
ALPHA  = 
BETA  = 


KY(  J  + 
,NE, 
KX(  J) 
KY(  J) 
»NF, 
KXC  J* 
RY(J* 
FXl-f 
FX2-F 

FYO/ 
*FY1 
FY2/ 

TER^l 
-(TEF 


GO  TO  4 


FX2)  GO  TO  b 
3) 
3) 

xo 
xo 

XI 

(D 
/( 

(D 

n* 

MO 


xoi*nxn2) 

DXni«DX12) 

xo?*nxi2) 

TEPM1*TERM2 
*(FX1*FX2)+TERM1*(FX0*FX2)+TERM2*(FX0*FX1)) 
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25 


5C 


100 


110 


200 


900 


D6R 

IF 

A  > 

P  ■ 

C  « 

IF 

IF 

DIF 

IF 

RfcT 

DIF 

PET 

DXO 

DXO 

DXl 

TER 

TER 

TER 

ALP 

BET 

A  s 

B  * 

C  s 

DIF 

IF 

IF 

ITI 

no 

IF 
IF 

ALP 
DIF 
PET 
DIF 
PET 
BET 
OIF 
RET 
FOR 
1  *. 
END 


IV  =  ? 

( A  d  s  ( r 

4,.AL 
4,«AL 
BETA* 
( AdS( A 
(  I  .LE 
uIS  ' 
(DIFDI 
URN 
DIS  « 
(jRN 

1  s  FY 

2  s  FY 
2  «  FY 
MO  e  r 
Ml  «  « 
M2  =  F 
HA  =  T 
A  3  »( 

4,«AL 
4,*AL 
BETA* 
DIS  = 
(DIFuI 
(DIFDI 
MFS  = 
TO  5 
(  ITEH 
(  I  ,LF 
IS 


, •ALPHA*DX*HfcTA 

ERIV)  ,nT,  1,  .AND,  ITIMES  ,EQ, 

PhiA«ALFHA 

PhA*bETA 

dETA^l, 

LPHA)  ,LT,  ,000001)  GO  TQ  10  0 

,  3  ,  ANP,  Y  ,tU,  0, )  GC  TO  25 

DSTINT(FX1)-UST1NT(FX0) 

S  ,LT,  n,  )  r;o  TO  20  0 

USTINT(FRX(3) ).nSTlNT(0,  ) 


1-^  YO 

2-FYO 

2-^  Yi 

X0/(DX01«D 

FX1/(UX01* 

X2/(DX02*D 

tRMO*TbRMl 

TEHM0*(FY1 

PHA*ALPHA 

PhA^dETA 

dETA*l, 

DSTINTCFYl 

S  ,LT,  P.) 

S  ,LT,  3.) 

2 


n  ^0  TO  50 


XU2) 

DX12) 

X12) 

♦TEHM2 

♦  KY2)*TFRMl*(FY0-»FY2)*rFRM2*(FY0*FYl)  ) 


)-DSTlNT(FYn) 
NO  TO  200 
RETURN 


,Lt,  1)  PRINT  900,  I 

,  3)  00  TO  110 

TOO  SMALL  TO  USE  PEGLLAR  FORMULA 

(FXl-Fxn  )*S(JRT(C) 


HA 
bis 

URN 

DIS 

URN 

A  a 

DIS 

URN 

MAT(«  IN  SUHKOUTINb  FRONDIS  ALPHA  /JAS  .ES?  THAN  ,00001  AT  POINT 

12) 


(F 


FRX(3)*SQRT(C) 

Yl-FYn)/(FXl-FXO) 
SOKT(RtTA*beTA*l,)*AfcS{FXl''FX0) 


SUBKOUTINE  FRONAM 

FHONAM  IS  TO  FINn  THE  POSITION  CF  THE  FkO\jT 

AT  THE  X  COORDINATE  LINES 
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10 


&c 


70 


100 


150 


175 


CUNHOiJ/A.l/AX,  AY 

C0MM0N/A?l/NI,NJ,MIl,,Nl2,NlvJ 
CUMMOiVA?iJ/NF,NF  1,NF?,NF3,  r'F'4 
COMMO.M/A30/ITtK,  ITFHT 

CrjMMO,>t/A40/Fr(?!5,6) 
CUMhON/A'iO/FRX(3n),KHY(^}n)  ,F  RS(3n) 

C0MMQiv|/A75/INu(2'^,6) 

CU^'MOl^l/A100/X(4),F  1(4) 

COMNOrj/r'Rl/LP 

CALL    F'ROripO 

DU    10  g     I=>5,rJIl 

A  I    =     (  I  •  3  )  *  A  X 

NOC    =    1 

NOP    =    2 

L    =     IiNPC  I,NUP) 

CALL    JFK(L,1JEHIV) 

IF    (AriS(ntRlV)     ,nb, 

LL    =    L+1 

CALL     UEl^(LL#nt:RIV) 

IF   (AdS(nePiv)    ,nT, 

Kb   =   L-1 

DO    in    J=l,4 

X(J)     s    FPX(KH+J^1) 

Fl(J)     =    FHY(Kri+J-l) 

CALL     INTFk(  AI,4,1,FF  (  I  ,M)C)  ,nLM,PUM) 

GO    Tn    70 

FF(  I  ,in|OC)     =    FRY(I.)  +  (FKY(L*1  ) -F  RV  (  l  )  )  *  (  A  I  -  "PX  (  L  )  )  /  (  FWV  ( I.- 

IF     (IThK    ,Lt,    1     .AMI).    LP    ,FU,     n)    PF-TIT    VOO.     l.L 


1,7)    CD    TC    !?0 


1,7)    GO    TC 


n 


1  )  -  F  -^  \  (  L  )  ) 


IF     (f'OC     ,UE 
MOP    =    Nur*i 
MUC    =     NOC+1 
GO    TO    5 
CUNTI  >IUfc 
DO    20  0     1=3,1 
IF     (irjD(I) 
MF     s     INiJ(  I  ) 
MFM     =     MF-1 


IND(  I)  )     yU    Ttl    inO 


III 

iLE, 


1)    ijO    TO    2L'0 


DO 
MM 

TF 
KK 
DO 
IF 


K=i,MFM 


17i 

=    K 

=    FF{ I ,K) 

=    K*i 

15U     HsKK.MF 

(  F  F  (  I  ,  M  )     ,  G  T 

=     1^ 
COMLJUfc 
TF     =     FF(  I  ,K) 
F  F  (  I  ,  K  )     =     F  F  (  1  , 
FF  (I,. 1M)     =     IF 


TF)     bfl     TO    1^0 


Mf) 
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IF     (ITFH     ,Lb,     1)     PPIM 
CONTINUb 


V^O.      I , (^ .Fr( I ,") , 


tl,VF) 


POO 

RfcTuRN 
900    F0RnAT(tJ<5H 
iTfc     ,I^,^2H 
950    FORmAT(^PH 

END 


WE  AHF  USln(i  LINEAR  I  MTF-h  pQL  AT  I  On  I  rj  FHPNAH  AT  C03-<l'iM. 

t^EGIN'  PJG  ViITH  POINT  ,  J?) 

AT  COORDINATE  LlNfi  .I^.^dH  Fr(,I?,4H)  =  .  F  1 1  .  t? )  ) 


SUBROUTINE  P03U.M  I  J.  IPOS) 
COMMOiJ/A?/AX,  AY 
C0MM0.J/A4C/^F(P5,6) 
COMMO  J/A75/INiJ(i?'^,h) 
Aj  =  ( J-1 )*mY 

IF  t  AJ  ,r;T,  FF  (  I  )  )  uu  Tu  1 
IPOS  s  0 
RETuR.Nl 
1  IN  =  IMJ( I ) 

GO  TO  (  in, 2u, 30,  40,^11  )  ,  IN 

10  IPOS  '    i 

RtTuR  J 
20  IPOS  =  lFTWfcN(Aj  ,LE,  FF(I.2),1,0) 

RETUP  J 
30  IF  (  AJ  ,Lt,  FF(  I  ,3n  GO  TU  20 

IPOb  =  1 

RETURfJ 
40     IF     ( AJ     ,LC,     FF( 1 ,4) )     GO    TU    30 

IPOS    =    0 

PETuP'J 
50  IF  (AJ  ,Lb,  FF( 1,5) )  bO  TQ  40 

IPOS  =  i 

RETUR'^ 

End 


SUBROOTIfiE    FRONPn 

C0MMn:j/A3/AX,  AY 

COMWON/A?i/NI  ,riJ,'JIl»Nl2.Nl3 

CUMhnN/A?2/NF,NFl,f;F^,IIKvi,'oh  4 

CoMM0.</A?3/aKT 

COMhON/A30/ITET, ITFKT 

CU^'H0N/A5U/FPX(30),FHY(>i0),HyS(3n) 

CUMMO  <i/A75/IND(2'^,ft) 

COMMOrJ/PRl/LP 
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INDW    =    U 

no   1    ICO   =   ;L,rjivi 
lNt)(lcn)    =   u 

DU    mo     if  1^    =    l,ijf--5 
DQ    mo     ICO    =    3,M2 
Ai    =     ( IC0-3)*AX 
Ul    =    FRX( IFK)-AI 
02    =    Ff<X(  irK*l)-AI 

Bi  =  Ji*n2 

IF    (Di    ,GT,     n,     .OH, 

1) 


suhkou 

TOONEK 
F  OK 

C  U  H  M  0  ^  J 

COMMON 
CUM^lOi^l 

CuMi-iniNi 

C  U  M  H  0  N 
CO^hOixl 

coHr,ni\i 

COMMON 


Tlr-:t     TOUMLfi  I  ,  J) 
KiNDb     Tub    VAU'ih 


Ul-     U,V,H 


THt     L'lFFEf'bNOt-     ^:UltATIO^S 
/A3/AX, AY 

/A'^/IJ(?4,27,3),  v(^'',;'7,3) 
/A6/Hl?4,t7,.3) 
/Al'i/K^^,-^?) 
/  A  ?  4  /  K  T 

/A3o/iTtP,  iTf-wr 

/A4U/^  F(^b,^) 
/AliiO/X(4)  ,(.-1(4) 

/  A 1  n  1  /  F  2  (  4  ) 


AT 
TU 


hF  i;sF 


Ton  ME  AW  TO  TUP  FH^MT 


J 
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'M  T, J,3)  .  I  .  J 


0)    FHiijT    son,    I, J 


fUl)    =   n, 
no  150    U»2#3 

IK    s    1-Ifv 
%{lii)     -     IK*  AX 
FKIN)     s    U{  l*If  ,  .J,3) 

r2<I')   =   s^(  l*iK,j,,i) 

150    r3(  I^  )    =    M(  UIK,  j,.l) 

CALL     IMPK(0,,'^,'',ll(I,.J,3),V(l,v.3),"(I,J.'^)) 


If     (Mn»  J.3)  ,LF.O.  )     CALL     I  M^F  (  U  ,  ,  ?  ,  1 ,  U  (  I  ,  J  , .-«  )  .  ;  M  ,  J  .  v5  )  ,  H  (  I  ,  J  ,  ?  )  ) 

RbTURN 
450    rORMAT(27h    h     IS    Ll-ISS    THAN    FtRC    AfJC    =     .fi^.'^.HM    AT    HOImT     (, 

1     12, H,,  12, IH)  ) 
SOO     FQRhAT  (lX,*t^F    Ubf^D    A    Rt-bUL^H    TOCNfcR    AT    pniu     (  *  ,  I  ^  ,  IH  ,  ,  I  ?  ,  IH  )  ) 
7Qn    FOKMATdJH    AT    POIiJT     (  ,  I  k  ,  1^  ,  ,  I  ?  ,  44H  )     IM    TDOilPR    Ln,J)     MS    GWEATbR 

ITHAN     9    AMU    =     ,12) 
7bO     FORMATdlH    AT    POInT     (  ,  U .  1^' ,  ,  I? ,  *  )        IN    TOD'^'f^P    l(l,J  +  l)     WAS    NJT    POS 

1ITIVC:*/100(1H*)  ) 
flOO     FQKhAr(it'*,72X,47Hl.'fc    UStl)    hQRIP'CNTAL     I  NTF -^PoL  AT  I  on     TO    FIMU    T^JJ'MtR) 
R50    FURMATdlH    AT    PUIiJT     (  ,  I  2  ,  1h  ,  ,  I  2  ,  Ih  )  ) 
900     FURHAT(Jf--H     IN    TUHNPH    UPPtH    POINT     IS    MlSSI>jn    AT     (  ,  I  ?  ,  1  H  ,  ,  I  2  ,  1  H  )  • 

1  74X,12(1M*)) 
95(1  FORNATCllH  AT  POINT  (,12,1^,, 12, 

1  *)   THt  POINT  (T,J-1)  IS  r'UT  MISSIKH  RUT  Is  TOUMbP*) 
960  FURNATdlH  AT  POINT  (,12, IK, ,12, 

1  *)   THt  POINT  n,J-l)  I^  MJT  MISSING*) 
END 


10 


20 


25 


30 


SUHROUTUS    CHbFH 

CHGFR    CHt-CKb     IF    FRX 

CUMHON/APO/TLeri 

CUMM0.VI/A22/NF  ,  Nf  1,  ^'|•  2,  Nl-  .S,  ^■F4 

CUMM0N/A?v5/AKT 

COMI'lOiM/A?0/KRX(30  ) 

CUMHOW/A=.2/MJ{:S0), 

IK    (FRX(^F?)     ,1-T. 

K    =    3 

F  RX(NF2)-TLFN 

FPY(NF2) 

FIjCKF  2) 

F  V  (  N  F  2  ) 


IS    LPSS    THAN    0    CR    ORP^Ttf'    THArg    20 


,FKY(3ll  )  ,FPS(3  0) 

F  V  ( .3 II  ) 

TLtN)     'iO    TO    4U 


TKFX 

TRFY 

TKFU 

TKF  V 

IF     (FRX(K) 

K    =    K  +  1 

GO    TO    ?U 

K^    s    NF2-K 

no    3n    j=i,KK 

JF    s     NF2''J 

rKX(JF-fi)     =     FKX(JF) 

Fr(Y(JF+l)     =     FkY(. iF) 

FU( JF+l)     =    FU( JF ) 

FV(JF*1)     =    FV(JF) 

FKX(K)     =     TPFX 

FRY(K)     =     TPFY 

FU(K)     =    TKFU 


TRFX)     (iO    TO    2! 
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iO  TO  100 


THFX)  CO  TO  50 


FV(K)  X  TRFV 

PHINT  200,  AKT,K 

GO  TO  10 
40  IF  (FRX(;^)  ,r,b,  '1.) 

TRFX  »  FRX(J)*TUf^N 

TRFY  «  FRY(3) 

TRFU  ■  Fl (3) 

TRFV  «  FV(3) 

K  «  NF2 
45  IF  (FHX(K)  ,Lt 

K  «  K«l 

GU  TO  4tj 
50  KK  s  K-1 

DO  60  Js3,KK 

FRX(J)  =  FRX(J*1) 

FRY(J)  =  FPY(J*1) 

Fo(J)  =  ru{J*i) 

60  FV(J>  =  FV(J*1) 

FrtX(K)  s  TPFX 

FRY(K)  =  TRFY 

FU(K)  =  TMFU 

FV(K)  =  TRFV 

PRINT  ?01,  AKT,K 

no  TO  40 
100  CALL  PEK?(1) 

FRS(2)  =  0, 

DO  150  1=3, NF3 

CALL  FRUNOIS(  1,0,  ,0,  ,  U  I  F ')  I  S  ) 
150  FRS(  1  )  =  FFiJ(  I-l)+DIh  Ulb 

CALL  PEK?{-1) 

RfcTUR.N 
?00    F0RMAT(29H    KRX    bXCFbUtU    Twhi>(TY 

i    18H    ANU    dECAMfl-    POINT     ,13) 
?Q1    F0RMAT(32H    FRX    WAS    LtrSS    THAN    ZFhO 
1    ISh    ANl)    dFCA,lF    PQINT     ,13) 

END 


AT     Tp<t     .F3,2, 


AT     TIMF     iff^t?, 


SUHhnUTlr.'t    KROMFH 

SUHROUTlr'E    FRONFH    IS    TO    F  p-iD    h    SUB    X    A^D    H    SUB    Y    AT    Thp    FrONT 

niMbNblOr'    MM(2),riIS(2) 

CUM|l0iJ/A3/AX,  AY 

C0MHnN/A6/h(?4,27,3) 

C0MMnN/A15/L(24,?7) 

COMMON/API/NI  ,NJ,fJll.Nl2,Nl3 

CUMMnN/A?i:/NF,NFl  ,NF2,NF3,Nh  4 
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T  ^F    r^n;^T    AT    priMT 


C0NMnN/A?4/KT 

COMhO.M/Avlu/I  Tt:W,  ITFRT 

C0MH0N/A4u/I-F(?b,6) 

C0MH0ivl/A5U/K«X  (3n  )  ,F  KY{^n  )  ,h  RS(jn  ) 

C0flNniJ/A51/l-HX(viP),hHY(«iU) 

COMMnN/AlOO/X(4) ,P  j  l^) 

COMMOW/PPl/LP 

KP  s  30D 

Wfc    ARt    THYINfi    TO    FIM'     THH     l''t[<  I VA  T I VF:  S    Of-"     H    AT    THF 

UE    FIKST    FIND    TUP    TA^'UbNT    TU    THh    FPC^T    AT    POINT     I 

Wfc    FIiMD    THIS    TANnfcNT    dY    USING    CIHIC     I  .MTtR^OL  A  T  I OM 

DO    300     Isi,NK2 

CALL    UFR( I^UhKlV) 

IF     (AdS([)fcPIV)     ,LT,     .OUODOn     EFrhlV 

SLOPF  =  -1,/nfcRlV 

SLOPF  IS  THt  SLOPfc  01-  THb  nUPf^AL  TO 

XX  =  SQKTd  ,*SL  1)^^*81  UPh) 

JPT  =  FRY( I )/AY*l, 

IF  (ABS(SLCPF)  ,LT,  3.)  HO  TO  lUf) 

FJP  =  JPT*AY 

nsT  -   Ays(xxw(Fjp-F(- Y(  1 )  )/sLnFP) 

IF  (KP  ,r-u,  0)  PPINT  9?a.  I#SLnPF,nST 

KH  =  IFThhN'(nsT  ,C,E,  ,?t>*AX,n,l) 

DO  99  Jlf'  =  l,? 

JpR  =  JPT  +  Ki)  +  Jl(>-1 

FJP  =  JPh»AY 

JPZ  =  JPP*1 

XS  =  FPX( I )+(h JP-FPY( I) )/SLUPE 

DIS(JINj)  =  AHS(XX*(FJP-KPY(  1  )  )/i)LCPE) 

XS     IS    THF     X    (:i;nKPir.iATt    of     THF    NjLPMAL    AT 

DIS  IS  THt  UlSTAMCF  HtT^Wf'  FRO^T  PCljT 

IPO  =  xs/Ax*:-;, 
IPI  =  iPn^i 

IF  (IPO  ,GT,  0 )  no  Tc  b 

PKINT  9d?,  I .  IPU, jPZ,SLOPEiFRX( J ),FRY(  I  ) 

CALL    i1YfcXlT(b) 


F  KONT 


Y  roOHnU'ATt  JPZ 

I  A;gn  THP  PJINT  (XS.J'^Z) 


b    IF  (IPO  ,LT,  Nl2)  no 


rn   lu 


fi)  bO  ru  IT' 


nUM  =  IFTHFN(JPU  ,Fu 

Kun  =  1 

RU  TO  4U 

10  IF  (L(iPn,jPZ)  ,';t. 

KST  r  I-? 

CALL  Fni5(  IPI,  Jp-^,-l.KSr,FX) 
FX  IS  THF  niSTAurt:  bbJHthN  IPI 
X(l)  =  IFi-FX/AX 

Fi(i)  =  n, 

DO  12    IK=^,i 

KN  =  IPO+IK-1 


Af-it  TM(^  F?J''T  CURVt  CRjSSINO  J-^? 
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x( Ik)  ■  KN 
12  F 1( IK )  *  H{KN, JP7,3) 

IF  (KP  ,FQ,  0)  PPINT  V40,  IPO, JHZ.rjD, VS, J1N,DIS( JIN) 

GO  TO  50 
15  IF  (L( lHO-1, JHZ)  ,GT,  0)  GO  TC  tO 

F IP  =  (  IPO-J)«AX 

CALL  .NFHPST(FIP,^JP,-1,^ST) 

CALL  FDIS( IPO, JP7,-li^ST,FX) 

X(l)  *  iPO-FX/AX 

FKi)  =  n, 

DO  17  IKs2,4 

KN  X  IPU*lK-2 

X  (  I  K  )  =  K  f  J 
17  FK  IK)  =  H(KN,  JPZ,3) 

NUH  =  4 

IF  (KP  ,GT,  0)  GO  TU  50 

IPOP  s  IPO-1 

PRINT  94n,  IPOP, JP7,^ JP»XSi JIN»CIS( jI  J) 

GO  TO  5  0 
20  IF  (L(IP1,JPZ)  ,nj,     0)     GO  TO  25 

KST  =  I-l 

CALL  FDISdPO,  JP7,l,KbT,FX) 

no  2?   iK3i,i: 

KN  s  IPU*lK-2 
X(IK)  =  KN 
22  FK IK)  s  H(KN,JH7,3) 
X(3)  =  X(2)*FX/AX 
Fl(3)  =  n, 
NUM  =  3 
IF  (KP  ,FQ,  0)  PPI^'T  94U,  IPl,  JF7,FJP,  VS,  Jl:>i.Dlo(  Jli>J) 

GO  TO  50 
25  IF  (L(  IPn*2,  JP7)  ,01,  0)  Gfi  TC  JO 
FIP  =  (  iri-«5)*AX 
CALL  .»)fcKFST(FIP,rjP,]  ,KST) 
CALL  FDIS( IPl, JP7,l,KbT,rx) 
DO  29  IKsl,i 
KN  s  IP04.IK-2 
X(  IK)  =  KN 

29  FK  IK  )  =  H(KN,  JP7,3) 
X{4)  =  X(3)*FX/AX 

fl( 4 )  =  n , 

UtjH     =  4 

IF  (KP  ,RT,  D)  LO    TO  50 

IPOP  s  IPO*<i 

PRINT  940,  IPUP,  JP7,F  JP,XS.  JIN,l.  ISC  jI  ^J) 

GO  TO  50 

30  KbO  s  1 
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100 


104 


122 


127 


y4i,    ^JP,Xa,  JIN, DISC  JI.nI) 


,Gf ,    rRY( r  )  )    UO    70    104 


NUM    =     4 

IF     (KP     ,FU,     0)     pt'INF 

40  DQ    41     IK  =  i,iMU(" 
KiN    s     IPJaNFU*IK-l 
X  (  I  K  )     =    K  N 

41  ri(IK)     =    H(KM,JF7,3) 
50     XA    =    XS/AX*i, 

IF  (KP  ,GT,  0)  bP  TU  V9 
PRINT  9on,  XA, ( IK,X( IK) ,  IK=l,NcM) 
PRINT  961,  (IK,F1(IK),  IK=1,NLM) 
99  CAL-L  INTFR(XA,NUM,1,HH(  JIN)  ,nLM,nijM) 
DX  =  XS-FRX( I ) 
DY  =  KJP-FRY( I ) 
RO  TO  POO 
IPT  =  FKX(  I  )/AX 
FIP  =  IPT*AX 

IS  =  -1 
IF  (FRY(I*1) 
IPT  =  IPT-»-l 
FIP  =  F  IP+1 

IS  =  i 

DST  =  AdS(XX*(FIP-FRX(I))) 

IF  (ITfiR  ,Lt,  1  .AND.  LP  ,ta 

KH  =  lFThtN(DST 

DO  199  JI,N  =  1,2 

IPH  =  IPT*KH*(JIM 

FIP  =  IPR*AY 

IPZ  =  IPR*3 

YS  =  FRY(  I  )*SLnpF-*(F  1P-FRX(  I)  ) 

DIS(JIN)  s  ARS(XX*(FIP-FKX( I ) ) ) 

JPO  =  YS/AY+l 

IF  (L(IPZ,JPO)  .nj,  0)  bO  TU  122 

KbO  =  1 

NUM  =  3 

IF  (KP  ,Eg,  0)  PPI^iT  V^lJ,  IPZ,  JPn.FIP,  YS,  JI 

r,U  TO  127 

IF  (L(iPZ,jpn-i)  ,c,^,    0)  un  tc  i^o 

KbU  =  2 

NU^^  =  4 

IF  (KP  ,GT.  n)  bn  TU  127 

JPOP  s  JPU-1 

PRINT  9';)0,  IPZ.  JPOP,  MP.  YS,  J  IN,  CIS  (  J  TM) 

x(i)   =  Fr(iPZ)>i 
Fi(i)   =  n, 

DO    12d     lK=2.NuH 
KN    =    JPU+IK-KHO 
X  (  I  K  )     =    K  N 
128    FK  IK)    =    H(  IPZ.Kfi,.!) 


0)     PRINT    923,     I,SUOPF,nST 


.LT, 
•1) 


,23*AX, TS,0) 
IS 


I ,  J  I  S  (  J  I  N  ) 
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130 


r,T,  0)  (iO  TC  133 


133 


137 
139 


140 
145 
150 

199 
200 


300 


0)  GO  TC  14U 

137 

IH, YS, JIN,CIS( 


IN) 


GO  TO  150 

If    {Hiy2,J^Q*l) 

KdO  =  ? 

NoM  =  2 

JPOP  *  JPO^l 

PRINT  950,  IPZ,  jPOP.t-  IH,  YS,  JIN,L1S(J!N) 

GO  TO  137 

IF  (L( IPZ.JPO*?)  ,nT. 

KBO  =  ? 

NUH  s  3 

IF  (KP  ,RT,  0)  GO  TO 

JPOP  3  JPO*ii 

PRINT  950,  IPZ,JPOP,^ 

DO  139  IK=1,NUM 

KN  =  JPU*lK-KHn 

X  (  I  K  )  s  K  N 

Fl( IK)  s  H( IPZ,KM,3) 

NUM  a  NUK*1 

X(NUH)  s  FF( IP7)*1, 

Fl(NUH)  s  0, 

GO  TO  150 

DO  145  lKsl,4 

KN  =  jP0*IK-2 

X  (  1  K  )  =  K  \ 

Fl(  IK)  =  H(  lPZ,Kr',3) 

rj  0  M  =  4 

IF  (KT  ,GT,  KP)  PfUNT 

XA  =  YS/AY*1 

IF  (KP  ,GT,  0)  Cn  TO  19V 

PRINT  960,  XA,  (  I  K  ,  X  (  I  K  )  ,  I  K  =1 ,  r^LM  ) 

PRINT  961,  (IK,F1(IK),  IKsl.NLM) 

CAUL  I^TEri{XA,MU^',l,H,1(  jlN),DLM,nuM) 

DX  s  FIP^FPXC  I  ) 

DY  =  YS-FKYtl) 

HM  IS  THF  VALUE 

=  JlSd) 

=  Jlb(2) 

=  rt2-Ml 

=  H2*HH(l)/(h1«nH)-Hl*HM(2)/(H2*rM) 

IS  M  SUB  N  AT  THE  KHUNT 

s    SIliN{l,,UX)*FM/XX 

s    SlliNd,  ,UY)*FN*Anb(SLOPF)/XX 

(KP   ,Ea,   0)   PPINT  y7o,   nx,i;Y,MK(i),Hn(  a)  .xx.Ffj 

=    FX 


951,     I-  IF,  YS,  JIN,niS(Jr>!) 


OF    H    AT    PuINTS    ARCVE     THE     --POMT 


HI 

H2 

DH 

Fri 

FN 

FX 

FY 

IF 

FmX( I ) 

FHY( I ) 

FHX     IS 

FHY     IS 

RfeTuRN 


s    FY 
THE    X 
THE    Y 


Dt-HIVATWH    OF    H    AT     ThF    F^^OigT 
DFHIVATIVE    OF     H    AT     ThF    FROivjT 
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92  0  F 

940  F 
1 

941  F 
1 

950  F 
1 

951  F 
1 

960  F 

961  F 

962  F 
970  F 

1 
980  F 

1 
985  F 

1 

2 
E 


ORMATC 
U  R  l-t  A  T  ( 
6H  XS 
QRMAT( 
6H  XS 
ORMAT( 
6H  YS 
ORhAK 
6H  YS 
URMAT( 
ORHAT( 
QRMArc 
ORMAK 
9H  Hc1( 
URMAT( 
IPO  = 
ORMA  r( 
7H  JPZ 
lOQClM 
ND 


4nx, 

8H  P 
=  ,F 
23H 
=  ,F 
6U  P 
-  ,F 
23H 

=  pf 
6H  X 
IHX, 
6H  F 
6H  D 
2)  = 
IX,* 
♦,  13 
i3H 

-  » 
*)  ) 


9HAT    C'OIi^T     ,  I2,10H       ilCPE 
On<jT     (,I2,iH,,I2,P>lH)     l^AS 


=   ,Fl: 
Missr 


12,6,5H    D1S(#  I1,4H)     =     ,F1?,^) 

NO    POINTS    Al-Kh    MISSING, 4X,7H    FjP     = 

12,6,5H    ill^i  ,  U,4H)     =    ,F1?,^) 

UIimT     (,  I2.1H,  ,  I2,2aH)     WAS    i^IS'^PIG 

l^.ft^-^H    DIS(  ,  I1,4H)     =     ,F1?,6) 

NO  POINTS  rtPKfc  f^ISSIi\r,,4X,7H  FJP  = 

12,6,?H  DIS(  ,  I1,4M)  =  ,Vl?  ,(-.) 


,7H  DST  =  ,F12,/) 
FJP  =  .F12,6, 

,F12,6, 

FIP  =  .Ft2.6, 


F1 2,0, 


,F12,^,4(4H 
FK,  12. 4H) 
H,6H  XX 
f ) ,  6  H  D  Y 


A  = 

4(4n 

N  =  ,F14, 

X  =  ,Fin, 


X(  ,  l?,4h)  = 
=  ,F1?.^) ) 
=  Fi2,6) 

=  ,M,D.6,9H 
»Fin,6,6X, 
GO  RbYONC 


■  1. 2  ,  6  )  ) 


.F10,6,6H  XX  s 
WF  ATTtMPThlJ  TO 
) 

IPO  WAS  Ll-bb  THAN  ZFKO  AT  Pnj 
12, '-^h  SLUKu  =  ,t-l?,6,7h  FPX  = 


H  -1  f  1 )  =  ,  F  J  0  ,  6  , 
6M  Fvj  =  ,F12,P) 

THE  iND  OF  THF  REnlON 


WITH 


T  ,12,1?H 
,Fl?,6,7H 


^ITH  ipn  =  , 12, 
FRY  =  ,F12,6/ 


SUBROUTU'E  l)EK(L,DFHI  V) 

CQMM0N/A24/KT 

COMMnN/A5  0/KRX(3n)  ,FKY(on)  ,F  PS(,in) 

DIMENSIU^^  D(?) 


XO  -    FRS(L-l) 
XI  s  FRS(L) 
X2  =  FRS{L*1) 
TtRMO  =  


X2  =  FRS{L*1) 

TtRMO  =  (Xl-X2)/( (XU-Al)*(Xa-X2)  ) 

TEHMl  =  (2,*Xl-xn-X2)/(  (Xl-Xn)*Ul-X?)  ^ 

TtRh2  =  (Xl-XO)/(  (X2-X())*(X2-Xl)  ) 

DU  3  ISLal,2 


IF  ( ISL  »uT,  1)  no  TO  2 
FXO  =  FKY(L-l) 

F  K  Y  (  L  ) 

FKY(L+1) 

3 

FRX(L-l) 

HKX(L) 

FRX(l.*l) 


FXl  = 
FX2  = 

GO  TO 

Fxn  = 

FXl  = 

FX2  = 
D(  ISL. ) 
DtRIV  = 
RETURN 
END 


=    FX0♦TEK^■0♦FXl♦TEK^'l  +  ^"Xi:*TFKr'2 
D(l)/C(2) 
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22 


24 


30 
31 


60 


Su 
r-b 
CO 
CO 
IF 
IF 
Mh 
WE 
Sri 

ro 

AS 
IF 

Srt 

MM 

CO 

IF 

Kb 

RE 

IF 

IF 

K 

K 

IF 

IF 

Kti 

R6 

K 

K 

IF 

IF 

KB 

Rt 

KB 

IF 

IF 

Kb 

Rk 

EN 


BRGJ 

PEST 

HHON 

HMON 

I  CI 

IDI 

«  3 

APE 

s  S 

5  M 

M  = 

(AS 
s  A 

=  M 

MIN 

(Fr^ 

=  I 

TURN 

ini 

(  iL) 

=  MM 
«  K* 
(K 
(FH 
=  K 
TURN 
=  MM 
s  K- 
(K 
(FK 
=  K 
T  U  R  iN 
s  I 
FRX 
FRX 
IS 
TURN 
D 


TINE  NERFSTi Al , AJ, IDIH.KR) 

FINDS  THE  NFA^tT  FHONT  PCINT  TO  TH=  COLUMN  Al 
/A?2/NF,NF1,NF?,NF3,NF4 
/A50/>'RX(3n),FHY(3O),FRS(Jn) 
R  IS  POSITIVF  THtN  FRONT  CURVE  IS  T^  T^E 
H  IS  NEGATIVF  THEN  F«OnT  CUkVE  IS  TO  THE 


KIGMT    or    HfcSH    i»OINT 
LEFT    OF    ;1FS>^    ■'OInT 


FINDING  THE  NEAHEST  POINT  ON  ThF  FRO:jT  TO  (Al.AI) 
0rtT((FKX(3)-.A|)«(FRX(3)-Al)*(FnY(3)-Aj)*(FOY(O)-AJ)) 
=4,NF2 

SQRT((FKX(f")-Al)*(FRX(M)-Al)*(FhiV(M),AJ)*(Fr<YC^)-AJ)) 
M  ,GE,  S"^)  GO  TU  b 
SM 

Ut 

V(M,1)  ,LT,  AJ)  GO  TO  22 

FTHEN(IDIf<     ,LE,     0,MM-1,KM-?) 

R  IS  1  WE  GO  T(i  THE  RIGI-T,  CThFRWISE  TO  THF  LETT 
IK  ,l.t,  n  )  GO  TO  30 


,oT,  NF2)  GO  TU  60 

V(K  )  ,LT  ,  AJ)  (lU  TO  24 

-2 


*1 
1 

.UT, 
Y(K  ) 
-1 


2)  r,c  TO 
,LT,  AJ) 


60 


TO  31 


FTl-EN(FHX(fM)  ,LE,  A  I  »  Mh- -1 ,  MM-?  ) 

(MM)  IS  LESS  THAN  AI,  wE  ISE  THE  pni'iTS   MM-1  ,  ^h  ,  mm^i  ,  >i  >i*2 
(MM)  IS  GRFATEP  THAN  AI  ^t     USE  THE  ^niMTS  MM-2 , MM-1 , MM i 1M+1 
THE  FIRST  FRONT  POINT  USFL  IN  THF  I  vjTERPOL  A  T I  ON 


Subroutine  nehfsty( i.aj, iuik.kh) 

common/a40/ff(?5,6) 

common/a75/inl(2'^,6) 

IF  lOIP  IS  POSITIVE 
IF  iniR  IS  NEGATIVE 
IF  (  I  ^rc  I  )  .GT  ,  1  )  G 


ThhN 

FWUNT  CUKVE  IS  A^^V"^  THF  MESH  POI^'T 

Then 

FRONT  CURVE  IS  r)  =  LOW  THF  MtSH  POINT 

II  TU 

in 
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MOC  =  1 

GO  TO  100 
10  IF"  (  IJIK  ,nT,  0) 

MQC  =  INDC I )*1 
15  MOC  =  NOC-1 

IF  (AJ  ,LT, 

no  TO  luo 

50  NOC  -    1 

55  NOG  =  NOCi-1 

IF  (AJ  .GE, 
100  KB  =  INU(  I.NOC*!) 

RbTURiM 

F:ND 


GO  TU  50 


FF  (  I  ,'400  )  UO  TO  15 


FF ( I, NOG) ) 
1 


^0    TO  55 


SUriK'OUTIht  PFKlOn(M) 

C:0MM0,>l/A5/U(2*;,i?7,3),V(24,?7,2) 

CUMrinr|/A6/H(?4,27,3) 

nOMHON/A?l/NI  ,MJ,Nll.rJl2»Nl3 
DO    1    J=1.NJ 


U(1,J,^^)    s 

V(1,J,M)  = 
H  (  1 ,  J  ,  M  )  s 
U(2,J,M)     = 

V  (  2  ,  J  ,  M  )  = 
H  (  2  ,  J  ,  M  )  s 
U  (  M  I  2  ,  J  »  M  ) 
V{M?,  J»M) 
H  (  N  I  2  ,  J ,  r ' ) 
U  (  N  I  3  ,  J  ,  t-i  ) 

V  (  N  I  3  ,  J  ,  M  ) 
H(riI3,  J,t) 
RbTURN 
f^ND 


U  ( I J I  ,  J  ,  M  ) 

V  (  N  I  ,  J ,  M  ) 
H(,jl  ,  J,f1) 
U  (  M  T  i  ,  J ,  h  ) 

V  (  N  T 1 ,  J  ,  H  ) 
H  ( l^n  1 ,  J  ,  M  ) 
=  U  (  3  ,  J  ,  H  ) 
=  V(3,J,H) 
=  H(3,J,H) 
=  U(4,J,h) 
s  V(4, J,H) 
=    H(4,J,h) 


SUBROUTlNb  PtH2(  IRFGiN) 

COMhON/A?0/TLtM 

COMMON/ A22 / NF,  NF  1,  MF  >^,MF 3,1 'F  4 

CUMhO.M/A5(j/FRX(3P),FHY(3n)  ,FRS(3n  ) 

COMMOi^/A?2/FU(30)  ,FVC3n) 

Pl:K2    fcVALUATFS    THE    FKONT    VALlitS     If^ 

IF     lE^tGIN     lb    0     WF:    DO    ALL    OF     PfeRi 

IF     IRtGIr.     IS    NPGATIVL    WE    [)')    SEClICN 


Af'CnRnAMCF"-     WITH    PFRIOPICirr 
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IF 
IF 


IBtGiN  IS 
( IdFGIN 


,GT. 


POSlTIVt  Wfc  DO  SbCTICN  2    ONLr 
P)  GU  TO  in 

Sf-CTION    1 


Fud)   » 

Fu<2)     = 
FU(NFo) 
ru(f.F4) 
FVd)    » 
FV(2)    « 
FV(NFJ) 
FV{NF4) 
FRS(l)     = 
FRS(2)     ' 
F«S(*.F4) 
FRX(l)     s 
FRY(l)     = 


FU(NFl) 
FU(NF2) 

*   ru(3) 

«    FU( 4) 

F  V(NFl) 

FVCNFiJ) 

=    FV(3> 

a  FV(4> 
FRS(3)*rRS(NFl)-FWS(NFJ) 
FPS(3)*rRS(riF2)-rRS(NF3) 
»  FPS(I.F3)*f  Hb(4)-FRS(«}) 
FRX(NFl)-TL.fcN 
FRYCNFl) 


IF  (IdBGiN  ,LT,  0)  RbTURN 

Si-CTION 
10     FRX(?)     s    FRX(NF2)-TLe-N 
F  R  Y  (  K  F  2  ) 
=    FRX{3)-»'TLbN 
s    FRY(3) 
=    hRA(4)4.Tl.hN 
a    FRY(4) 


FRX(?)  s 
FRY(?)  = 
FRX(NF3) 
FRY{\F3) 
FRX(rJF4) 
FRY{MF4) 

END 


SUBROUTlr.b  INT&h(XX,f!UNPiNUMPEP,  AINI.AIN?,  Ap'3) 

INTER  PtPFOKMS  LAGRANGE  I  NTERF OL A T I C N 

CO^^Mn.N/A100/X(4)  ,F1(4) 

COMHON/A101/Fi:(4) 

CC^*H0N/AliJ2/F3(4) 

NUMP  IS  THE  NUMbER  0^  I NTEWPOL  A  T  I  KG  POINTS  U'^EP 

MUMbEH  IS  THE  Mf'HFP  UF  FUMCTIONS  DESIRED 

REAL  NUM 

AINl  s  0, 

AIN2  =  0, 


REAL 
AINl 

AIN2 
AIN3 

DO  IC   ^ 
HUM  =1,0 
=  1,P 


NUM 
=  0, 

0 
=  0. 
KL«1,NUMP 

^  <   n 

DEN  =  1,P 
nu  4  JL=1,NUMP 
(KL  .FCJ.  JL) 
■X( 


IF  (KL  ,FCJ, 

NUM  =  MUM* (XX    

DfcN  =  Dbri*(X(KL)-X( 


no  It'  4 

JL)  ) 

JL)  ) 
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CONTINUE 

IF    (AdS(ntN) 


PHINT    9l?n,DtN,NUfi,KL 
GO    TO    60 
DNUh    s    HUH/Df:N 
AINl    s    AIN1*F1(KI,  )*U!MUM 
IF     (NUMbFR-ii)    1U,8,9 
8    AIN2    =    AIf'J?-^F2(KI.)*;jf'Uf- 
no    TO    lU 
AIN2    =    AIN2*F;c;(Kl   )*l"jMuM 


r^t,    .Duununi)  r,c  lo  7 


PRINT  99B 

CALL  tXiT 

950  FURhAT(0H  

960  FuRMAT(lX,4(4h 

961  FuRMAT(1X,4 (4H 

962  FURMAT(1X,4(4H 
998  F0RhAT(5n(lH*) ) 

r-ND 


F1  (,  I^,4H) 
F  ?(  ,  U,4ri) 


.Ms 
=  ,  F  1 4  ,  fc  )  ) 
=  ,F14,fa)  ) 
=  ,1-  14,b)  ) 


,F1P,5,7H   TIMb  ,12) 


SUHROUTH'b  HFLAbLE 

SUHRnuTirt  HFLAm.^  RhuIbTRlbHTFb 

rnM,-in,M  /  Aon  /  Ti  t^M 


COMi-10iM/A?0/TLEN 
C0MHnN/A?2/NF,MFl,rJF2,UF>i,r'F4 
COMMON/ A5  U/FRX(30  )  ,h  KYCOiH)  ,h  PS  ( 
C0MM0N/A32/I-  U(,"^0  )  ,FV(.VJ  ) 
CUMMON/AlOn/X(4) ,F1(^) 
C0MM0In|/A1D1/F2(4) 
CUMI-iniJ/PPi  /LP 


THE  POINTS  ALONG  THF  FRONT 


iO) 


COii|"i>jH/  n  X  ij  A.  f  r 

CUMMON/PPl/LP 


U)  ,TFRb(3n),QFR?(-^0) 
(TK V(l)  ,TFRY(t) ) 


i^Ui  ii'iui>i/  rf^  J./  ur 

niML-NSJlUrj    Kb  (30) 

niMfc^'SIUf    TFRX(3P),Tt- RY{v5u)  #T 

DIMLNblUrj    TMJ(3(j),TFV(^0) 

FgUIVALbMCE  (TFU(1),TFRA(1)), 

HbLS  =  (FRS(NF3)-FRS(3) )/NF 

PRINT  500,  DELS 

TFRX(3)  =  0, 

TFPo(s5)  =  n, 

DISTAMCb  IS  MbASliRFD  KKUM  ThR  PCIKT  (j,Y(J)) 

DO  1  K=l,3 

X(K)  =  FPS(K-fl) 
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9 
10 

15 

20 


25 

30 


40 


a  FriY(K*i  ) 

NTfR(0,,3,l,TPHY(  J),i)UH,DlM) 
E  JUST  POur-U  TI-KY(4) 
a4,NF2 


TKRb(  I-l)*l'cLS 
;F2 


NF35 
•  LT, 


no  TO  b 
TFRSC  I  )  ) 


GO  TC  d 


FKK)  a 
CALL  I 

we  hAv 

DO  2  1*4, NF2 

TFWS(l)  = 

KB(J)  =  1 

ro    5  1=4. 

K  s  1-4 

K  «  K*l 

IF  (K  ,uF, 

IK  (FHSlK) 

Kb(  I  )  =  K-2 

DO  10  1=4, NK2 

DO  9  K=l,4 

L  ■  Kti(  I  )*K-1 

X(K)  s  FRS(L) 

Fl(K)  s  FrtXlL) 

F2(K)  =  FRY(L) 

CALL  INTfeH(TFRS(l),4,2#TFHXlI),TFhY(I),DUvi) 

^.b  HAVE  JUST  FnUK'D  TKKX  ANU  TFRY 

DO  15  I=1,N'F4 

OFRS( I  )  =  FKS(  I  ) 

DO  20  1=3, NF2 

FRS( I )  = 


FRXi  I )  = 
FHY{  I  )  = 
FRX(NK3) 

FHX(r;F4) 

FRY(r.'F3) 
FRY(Nr4) 


TFKS(  I) 
TFHX(  I  ) 
TFHY(  I  ) 
=  TLEN 

=  FRX(4)*TLt:H 
s  FRY(3) 
FRYM) 


30     l=il,N'F2 
25    Kxi,4 
B     Kd(  I  )*K-1 
OFKSIL) 
»    FU(L) 

s   rv(L) 


RNFHS 


DO 

DO 

L 

X(K) 

FKK) 

F2(K) 

CALL 

Wt     H 

DO    4P 

rod) 
FV(  I ) 

CALL 


|NTEK(FRb(I),4,i;,TFU(l),TFV(I),njM) 
AVE    JUST    FnuriQ    Fli    AND    FV 


i=:^,rjF  2 
=    T  F  u  (  I ) 

=     TFV(  1  ) 
PEK?(0) 
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RtTURiJ 
500    FORMAT(lX*SUHKOUTirjt 


HhLAfcil.  h 


DHLS    =    ♦,F1P>,7) 


FO 
1  7H 
930  FORHA 
950  FORMA 
FURhA 
END 


SUBKOUTIME  MYPRr.T2(H) 

C0MI-10iN/A5/U(?4,;:;7,3),  V(^4,;^7,3) 

C0MMON/A6/H(24,je7,3) 

COMMOiM/APl/NI  ,NJ,NIl,NliiilMl3 

PRINT  460 

PRINT  4dn,  (J,  (M( I, J,M),  1=3, N 

PRINT  461 

PRINT  4dr,  (J,  (li(  I,  J,N)  ,  1=3, N' 


PRINT 
PRINT 
PRINT 
PRINT 
PRINT 
PRINT 
PfcTUPi^l 


460 

4dn, 

461 

4dr, 

46? 

4dn, 


(J, 


1=3, Nil), 
I =3, Nil), 

I  N  U  )  , 


J  =  1,'!J) 


1=3, 


)sl,'lj) 

Jsl,Nj) 


PtTUPi^l 
460  FOPMAT(//60X,lHh/3H 
1  1H6,11X,1H9,11X,2H1 


(  V  (  I  ,  J ,  h  )  , 

I  =  ,«3X,1H3,llX,lh4,llX,lH5,llX,lh6,llX,lH/,llX, 
(|.10X,?H11,10X,2H1?//) 
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461  rOPMAT(//60X,lHL/3H  la,&X,lh3,HX,lh4,liy,lM''i11V,lH(S,ny,lH7,llX, 
1  lH8,llX,lhV,llX,2HlU,10X,?Hll,inx,2Hl?//) 

462  rORMAT(//60X,lHV/3H  I«,6X,lh3,llX,lh4,llX,lH5,llX,lH6,llX,lH7,llX, 
1  lHb,llX,lM9,llX,2H10,10X,2Hll,10x,2Wl?//) 

480  rORNAT(iy,I2,in(Fll,4,lx)/3x,10lF.ll,4,lX)) 
997  FOWMATdH  ) 
END 


100 


SUBRO 

COHMO 

COMMO 

COMHO 

COMMO 

CQMMO 

COMMO 

YMIN 

SIZE 

TLEN 

SFX  I 

DISX 

DISX 

YMiN 

S^X  s 

SFY  = 

LF  « 

DIVX 

DIVY 

DIV  = 

CALL 

CALL 

CALL 

CALL 

AHR  = 

CALL 

CALL 

PRINT 

CALL 

RETUR 

FQRMA 

END 


UTINE 
\/A?0 
•'J/A21 
N/A22 
i>l/A23 
N/A9Q 
W/PL/ 
IS  TH 
IS  TH 
IS  TH 
S  TMfc 
IS  TH 
=  SIZ 
«  6. 

TLEN 

(TLE 
N  F  ♦  1 
«  25, 
a  25, 

10*K 
AXIS( 
AXIS( 
LINE( 
SYMHO 

AKT/ 
NUMHfc 
PLUT( 

100, 
CONT( 
N 
T(VH 


MYPLOT 
/TLEN 

/Nl  ,'^J,NlliNl2,Nl.3 
/NF,\F1,NF2,NF3,MF4 
/AKT 

/l-RX(30),FHY(30),hRS(30) 
SlZbX,SIZbY,TLfcNY 
E  MINIMUM  OF  Y  TO  RE  FLCTTEIJ 
F  LENGTH  OF  THE  AXIS  IN  INCHES 
E  NUMBER  Of-  COORUlNATt  LINES  OM  THE 

SCALE  FACTOR  FOR  ThE  X  AXIS 
F  LENGTH  bbTWEEN  PLOTS  IN  INCHES 
FX*2,5 

/bIZEX 
NY-YMIN)/SIZEY 


AXES 


4 

4 

ANRE/( 
0,.0,  . 
Oi.O., 
FRX( 3) 
L(.75, 
60  , 

P(,75, 
DISX,C 

AKT 
NI.NJ) 


(LENGTM*NUMHtR    O^     CnCRJIi>lArP    Ll^tS    PER 

,iHx,-i,sizFx,o,,n,,sFx,uNy) 

,1HY,1,SIZEY,Q0,,YMN,SFY,DIVY) 

),FRY(3),LF,i,i,ii,n,  ,sFx,  y"rj»'=;FY) 

,■=-,0,  ,28.bHHOufiS,0,  ,b) 

,5.5, ,28, AHR, 0,,^) 
C.  .-3) 


TIT) 


AT  TIME  ,F12,7,*  WE  CLMPLETED  A  "L^T*) 


SUBROUTIf.E  MYbXlT(N) 
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50 

100 
200 


C0MMni\j/A25/IPhIM,IPL(jT 

COMMOi\l/PF<l/LP 

IF  (LP  ,1-U,  0) 

IF  (LP  ,EU,  0) 

IF  (IPLOT  ,Le, 

CALL  PLUT(5,,Q 

CALL  PLUT(0,,C 

PRINT  200,  N 

STOP 

FaRHAT(19H  Wb  AKf"  NQw 

F0RhAT(6H  STOP  ,12) 

END 


Pf'U'T 
CALL  I 
0)  r,o 
,  ,-3) 
,»«99) 


10  0 

iYP|^NTl(3) 
TO  50 


fcXITiNG) 


SUt^KOU 
CUMMO.^J 
CUMMOU 
CQMMO.M 
DIMtNS 

EQUIVA 
EXTtRN 
DU  1  I 
DO  1  J 
HT( J,  I 
KUIH  = 
NNI  = 
MNJ  = 


E 

/H 

2/ 

0/ 

H 

I 

CF 

FX 

NI 

NJ 

M 


LOtxT  (ri  .hJ) 
(24,27,>1) 

NF,NF  l,MFi',Nf-  >i,NF  4 

i-HXCiO)  ,FKY(v5n)  ,KPS(3n) 

T(27,?1),(;l(5) 

TITLt(2)»LAHhLX(2),LAHeLY(2) 
(HT(1 ,l),h(l,l,2) ) 


(  I  *  2  ,  J ,  3  ) 


XA 

xy 

YA 

YB 
XG 

YG 
NCL  = 
DO  2  I 
CL(  I) 
ITITLt 
ITITLt 
LABbLX 
LABtLX 
LABLLY 
LABELY 
KP  =  Q 
CALL  C 
.  ITITU 


TIM 
/A6 
/A2 
/A^ 

ION 

lar.i 
LErj 

AL 
=  1. 
=  1, 
)  = 

27 

NI 
NJ 

I 

I-l 
I 

J-1 
I 

4 

=1,NCL 

=  ,l'^i^52*I 

(1)  =  lOHCONTOUi^S 

(2)  =  n 

(1)  =  6hX  AXIS 

(2)  =  0 

(i)  =  6hY  AXIS 

(2)  =  n 


ONTOUK(HT,KUIM,MNJ,NNI,f'Nw,KNI,XA,X-j,Y^,Yd,Xn,YG,NnL,CU» 
E,LABbL.X,LARFLY,FX,FX,KP) 
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YSf  =  (JJ-D/YG 

CALL  Li'Ne(rKx(i),rPY(3),Lr,i,o,i,o.,xsr,o.,Ysr) 

CALL  PLOT(n,,-ll. ,-6) 
CALL  PLUT(10,5,1. ,-3) 
RfcTuRN 
END 


FUNCTION 
rx  e  X 

RfcTURN 
END 


FX(X) 


SUHROUT 
1  M,N,MM 
DATA  TW 
COMMON/ 
COMMON/ 
COMMON/ 
lHX»hY,X 
COMMON/ 
COMMON/ 
DIMENSI 
DIMENlil 
Z(I,J) 
I  VARIE 
J  VARIE 
Z  HAS  D 
MXN  IS 
MMXNN  I 
XA,XP, Y 
OF 
XG  IS  T 
Yli  IS  T 
NCL  lb 
CL(  I)  A 
ITITLt 
LABcLX 
LAbELY 
THE  X(  1 
LIKbWiS 


iNt  COnT 

,  N  N ,  X  A  ,  X 
0PI/b,2B 
HOLAHC/H 
INDICES/ 
XYbNUS/X 
S,XSS,YS 
CLEVELS/ 

CAviN/in 

ON  Zd) 
Of"'  CL(l) 
IS  THE  0 
S  BETWEE 
S  oETWbE 
K^ENSION 
THE  SIZE 
3  THt  SI 
A,YH  ANE 
X  AND  Y 
HE  WIDTH 
HL  HbIGH 
THE  NO, 
WF  THE  C 
Ctl.'JTA  KmS 
CONTAINS 
CONTAINS 
)  ARE  AS 
b,  THE  Y 


Ol)R(Z.^DIM, 

fa,YA,Yb,xn,YG,NCL.CL,lTITLE.LAdFLX,LAR£LY.FX,fY,AP) 

3ia53l'7l7yt>fl/ 

s.Rn, Tms,tho 

MPOW,NCOL.MMRO^,NNcnL,KPOL 
MlN,XNAX,YhIN,YMAX,XSTZE,Y3lZF» 
, YSS,KXA,F  YA 
NLVLS,NI.V,CLbVEL('50) 
IM,nUM(4u35) 


RTINATt  AT  POINT  X( J),  Y(  I) 
N  1  AND  M 
h    1  AMU  N 

(KniM, , ) 

or    THb    CALCULATkD    X-Y    GRID 
ZP    OF     THt    bXPANDfcn(PY     IMTF^DQI  ATION) 
THE    MINIMUM    AND    MAXIMUM    VM.uFS 


X^Y    HRID 


OF    TMb    li^^APH       IN     INCHES, 

T    OF     Trie    liRAPH     IN     IfJC'ES, 

OF  CONTOUR  LEVELS 

ONTOUH  LEVELS 
THE  HLUT  TITLES      IT  SHOjLD  ENJ 
THE  LAhbLLiNG  ALONG  THE  X  AXIS 
THE  LAHELLING  ALONG  THE  Y  AXIS 

Si'MED  TO  HE  bOLALLY  SPACEH,  aND 

H), 


WITH  ZERO  WJ"<D 
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10 


FUMCTini.j  Til  f  fc  PI  OTTHL  ALC'i 
FUMCTirtfM  Tn  l^fc  PluTTh-i.  ALC' 
CARTFSIAM  CODrO)  p  ATb  S  Af-F  !i 
PULAP    roORDlNATcS    ARF    llSFC, 


0)    an    TO    i 


FX     15    THF 
FY     IS     TH[ 

IF  KP  =   n 

OTHt^'Jlbt 

IDlN  =  KDil 

MKOWsi-l 

NCOL    =     'M 

MriRUW     =    t    1 

NNCUL    =    f  iJ 

KPOL  =  0 

IF    (KP    ,rt 

XHlfM  =  XA 

XHAX    =    Xt' 

YMli-J    s     YA 

YMAX    s     YF 

GO    TO    ? 

XMIN=-Yd 

XMAX=Yri 

YMllN  =  XMir. 

YMAX=XMAX 

THO=XA 

PU  =  YA 

THS=(  XH-XM)/f  L')aT(M(vi;uI 

Pb=(Yd-YA)/H  UATdlMi^'i^. 

KP0U  =  1 

CUNTI.>Jllt 

XSIZEsXG 

YSIZfc     =     Yu 

NLVLS    =     lAf-'bCi^TL) 

CALL    PLUT(0, ,-,3*(ll. 


J    T^iE 

K- 

■AXIS, 

'J    THE 

Y-AXIS, 

i^D 

-1 ) 
1) 


Yoi£f  ),;i) 


^.f- 


u)   f',0   ro  I. 


IF     (NCL 

CLfcVFL    =    , 

HX    =    Z 

L    =    0 

DU    d     1=1. 

nu    7    J=l, 

L    =    L  +  1 

IF     (Z(L) 

CLfcVFL    = 

IF     (  Z  (  L  ) 

HX    =    Z(L) 

CONTIi^Llb 

LsL-t'+I  JIM 

HX  =  (HX-CLtVtL)/FLOA  T(i.L.vLi-l) 

DU  in  1=?, NLVLS 

CLbVfrLC  I)=CLEVhL(  I-U+HX 

CU^T  I"^Uc 
nu    TO    ?U 


■  ICUL 

DPUW 

•  GL, 

C  L  r-  V  F  L  ) 

/(L) 

,l-c, 

HX)    no 

i,o    Tn    5 


TH    7 
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12 
lb 
20 


30 


100 


200 


PO  15  1 

CLtVFLC 
HXa{XMA 

HY3 ( YMA 
XS«( XMA 
YS» ( YMA 

rxA=rx( 

FYA  3  F 
XbSsXU/ 
YSS=YU/ 
CALL  IN 
DU  30  N 
CALL  SC 
CALL  LA 
IF  (KPP 
XCfcNTsX 
YCtUTsY 
PZbHO=X 
RMAXrXC 
THtTAls 
IF  (THt 
IF  («Zb 
CALL  AH 
THH  =  T 
K  =  ? 
IF (AB3( 
DO  100 
SIMHsS 
CUSThsC 
XO=XCtN 
YO=YCtM 
X 1  =  X  C  c  N 
Yl=YCtN 
CALL  JA 
THHsTHt- 
CJNTI'U) 
CALL  JA 
CALL  JA 
CALL  JA 
CALL  JA 
CONTPMU 
PbTURN 
END 


=  1 
1) 
X» 

X- 

X- 
X- 
Xf 
Y( 

{y 

(F 
TE 

LV 
Af- 
dF 

L 

SI 

bl 

CE- 

XH 
FA 
KP 
C( 
HO 


.NL 

AM 
YMi 
XKI 
YMl 

IN) 
YM  I 
X(X 
Y(  Y 
KP( 
=  1» 
(Z, 
L(  1 
.EJ 
U* 
ZF* 
igT* 
T 


VLS 

CL(  1  ) 

N)/Fl  LAT(.iCJL-1  ) 

^i)/^  LuAT(/1NUri-l  ) 

N)/!'  LOAT(NtlL(ii_-l  ) 

N)/FLnAT (MMhDrt-l  ) 


N) 

MAX 
MAX 

Z,I 

NLV 
FX. 
TIT 

.  0 

it? 
.5 

Hu/ 


)  -  F  X  A  ) 

)-FYA) 

Dill) 

L*^ 

FY) 

LF,LAHbLX,L/«bFLY,Fx.FY) 

)  (in  Tu  ^'JO 


YP 


1  ,GT,  TwnPI)  THETAls  T'v(.,PI*A"'i)n(  TmFTA.TwOPI  ) 

.LiT,    1.)    CALL    ARCCXCFNT,  YCF:jT,r?ZFh(n,THn,THFTil,l) 
XCfcNT,YrfcrJT,Kl-1AX,Thn,ThFTAl,l) 


TMbTAl-(  TlJuFl*THL  )  )  ,LT  ,  ,  (Jfil  )     w     =    1 

1=1, r^ 

ir  (THH) 

U?(Tr<H) 

T-..KZtRU*CnsTH 

T*r<ZtRO*olrjTH 

Tfr(r-AX*rUf^TH 

T  +  KVAX^SriTH 

bHLlM(Xn,YO,Xt,Yi, ,1) 

TAl 

SMLlh(XCu'iT,  YUhi^iTiG  ,  ,YCFi\T,  ,2'^) 
bHLlN(XCt"T,Y(:tMT.Xi:tMT,  YSr/b,  ,?5) 
iiHulN(XCt*  T,  Yi:tMT,XCtMT,0  ,  ,  ,2'^) 
3i-Lli'(XCb"T,  Y(:e.M,XSlZE,  YCLMT,  ,?5) 


SUHHOJTlft    LAoFL(lTllLL,l.Ai'tlX,LAi,FLY,Fx."Y) 
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10 

11 


12 


14 
15 


C 

C 

iX 

c 

c 

D 
D 
J 
X 
P 
G 
D 
X 
I 
X 
X 
G 
I 
J 
Z 
U 
C 

c 
I 

p 

X 

w 
e 
w 
c 
w 

D 

c 
J 

Y 
P 
G 
D 
Y 
I 
Y 
Y 
G 
I 
J 
Z 
1^ 

c 


UMNO  4/ 

b.XSS, 

OMMOivl/ 

IMfcNiiil 
ATA     IS 

=  0 

=  XA 

=  0 

=  XSI 
0  10  I 
G  =  XS 
F  (Xij 

=  Xd 
G  =  XS 
U  TO  5 
}■     (XU 

=  J*l 
(J)  = 
fc  DPAW 
ALL  SY 
ALL  -JU 
^  (X  , 

=   y  ij  + 
=   x*x 

h    APt 
ALL    n 
t    Ar^t 
ALL    TI 
t    APE 
U    12     I 
ALL    SY 
sO 
=    YA 

=   n 

=  YSI 
Q  iiO  .  I 
G  s  YS 
F     (Yb 

=    Yri 
G    =    YS 
Q    TP    j 
F     (YG 

=    J  +  1 
(J)     = 
t     A(?k: 
ALL    iJY 


TfriiP/Zdul  ) 

XYMhUS/XA.Xf  #  Y«,  Yr,  >bI2F-*  YSIZ^.Hx,  HY, 

YS,YSS,rXA,fYA 

iNJICf^b/t-  ,.N,.1",l''i^.K'-'UL 

CLtVtL3/i\n|_,  ,mL  V»CL('>d  ) 

OM     I  r  ITLt(l)  ,LAi.iLl  X(l)  ,LA!U  LY  U) 

/18/ 


Zt-,9 

=  1 ,  M IM 

b*(FX(X)-FXA) 
.LT,     H)     Ln    Tu    •+ 


iZt 


,LT,    P)    (.r'    T.J    1(1 


xr 

A  h!  p  U  W 
HPUL (X 
llbLRCX 
L  r  ,  X  H 
.P 


b    Tn;iPTt'tiv    wjTh    ^U^H,^hS    aT    Thc 
bj-,07,,ii'.lSiU,,-j.) 

I-  - ,  -5  ? , ' .  <♦  (' , ,  X  "^  f  y ,  0 ,  ? ) 

)     uf    TO    il 


■Oirnn     UP     THP;     r.RAr'H 


ru.     TML     X     MX  Is 

, -.(^'L),  AMZt.  ,14,0,,I.AMtLX) 

I.    A    TITIb    Al     IhP    TCP    OF     THE    3"A'''i 

#YSTZr+,l'?.XbIZF,.iil,0,,ITIT,r) 

u    AFl^nivS    Al     Tl't     TOP     (-F     TrlP    '•,^up\^. 


LAbPLl- 
rLt(U, 
PliTT  IM 
TI  L(U, 
PL^TTIN 

=  3  ,J 

ilPOL(Z(  F  )  ,  YSi/L-^.ny,  ,12,  IS,!^!^.  ,-1. ) 


ZP-,'.d 

=  1  ,  M  i-l 

b*(F  Y(  r  )-FYA J 

,LT,     lO     uT:    TU    1^ 

lit 

,L  r,    P)    un    TU    «^U 

Yr 

UPaWH-iG    APK'flWS     Tui-pTHF-h     WITH 
|1PUL(-,1  4,  YH,  ,<:;b,  l,S,i^70,  ,-1) 


M.jMdP^'-,    AT    Tt't    L^'PT    nr    ij-<HHh 
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CALL  NL'MRtB(-,S&,Yn-»,l^,,l«iY,n,;^) 

IF  (Y  ,tr.,    YB)  on  TO  21 

P  =  Yi4*  ,9 

20  Y  =  Y*YS 

Ut    ARfc  LAdELLiriC,  THK  Y  AXIS 

21  CALL  TITLt(-,d,U,  ,  YSIZH,  .I'^.'JU,  ,L  ARtl  Y) 

DO  2?    Isl»J 

WE    ARt    PI  AflNi,    ArHnWb    t1i«.     THb    RIGHT     S  IJF    0."     T^^h    JI^APH 

22  CALL    SV;-U^0L(XbI2P*,U/,Z(  I  ),  ,75,  I-^.tJO,  ,-I) 
YI     =     ,5*vS!Ze*,l*KL0AT(h<:L) 

DY    s     ,4 

WE    ARE    LAdELLTHj    T^h    COOiauUW    CLt'VFS    ON    THF    pIOHT    OF    THE    GRa^h 
CALL    .^UMHEF(XSTZf-*,i:»y  I  .  .  ii^^  i  MCL  *  0  ,  ?h  I  2  ) 
CALL    SYMROUXblZk-*,  Jll,Yl  ,  .>'d,14hCCNT'^UR    L=VtLo,0.14) 
Y  I  s  Y  I  -  r)  Y 
DO    24     I=l,rjCL 

CALL     SYMf'0LUSIZr:*,o,Yl*,(lb,n,2u.  l-l,Ji-1.  ) 
CALL    .■iUHDEP(XblZP*l,",  Yl,0.205,LL(  I  )  ,  0  i  "s'^  ^1  '3  . '? ) 
24    YIsYl-DY 
RcTuRiN 
END 


T      =     HA+AA 

UU    13    J  =  ?,M 

CALL    KIT( J,T,nX, AC, 

12  AM(  I  ,K)=7Fl.i>.(  aY) 

K     s     K*l 

XY     =     XY*XS 

IF     (XY     ,Lc,     T)     b*^    T 

13  TsT*i-'X 


1,A2) 


T)     b*^    TO    i^' 
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14 


15 
16 
17 


18 


19 


20 
21 


2b 
30 

50 


IF  (K  ,UT,  imM)  bO  TO  1^7 

AM(I,.<)  =  ZI-U>;(XY) 

K  =  Kfl 

XY  =  XY+XS 

fiO  TO  14 

CONT  I.NUt 

IF  (M-Mrl)  l7,Jil,'^o 

no  25  1=1,  Mm 

DO  18  J=1,M 

Z{ J)  =  AM( J,  I  ) 

CUNTIiNlJt 

K    =    1 

XY    =     YA 

T    =    HY*YA 

DO    20    J  =  ?,M 

CALL    KIT(J,r,hY,AO,Al,A£) 

AM(K,  I  JsZHUNCXY) 

K    =    K  +  1 

XY     =     XY-fYS 

IF     (XY     ,1.  b,     T)     UP    TG    iV 

TsT*HY 

MM) 


IF     (K     ,uT, 
AM    (H,!)     = 

K     =    K  +  1 

XY    =    XY*YS 

RO    TO    ?1 

CONTI-lUb 

HkTUP.J 

STUF12 

END 


i,p 


TC. 
(XY) 


^.i> 


SUHhOUTlr't:    brA"(  ri,Fx,H  r  ) 
AM     IS    THf-    MATHIX     T'l    f'L-    U'),,  Tijnwi,  J ,     Ml     atJC 
CL(NLV)     I'i     Tl't.    C();jHj|iK^    Lt-Vl-L, 
THb    M        (X,Y)        VALUt-b    ijf"     D.jfc    Cui'TOLR    LINE 
THtY     AHF    AVAILAHL^-, 

DIMtnSIUf      A.i(i) 

COMHO:-l/CLtVhLo/;-i"Li  ''l-^*  i-L  C'-J  ) 

CoMMO,\i/ir.iUIUtS/DMM(c)  ,NT,NT,KFnL 

CuMMON/CAVlN/uIi!,  I  X  ,  1  <  .  1  HX  »  I  D  Y  ,  I  SS  ,  KP  , 'i ,  H  ;  ,  J 
1  I  ^JX(M)  ,  iriY(c)  ,l?t(.  ('SOn  )  .X(I6o;<)  ,  Y  C  160")) 

TYPt     lUTI-UFK    Kf-CJl." 

DATA(  IMY=U,i.i,i,n,-l  ,-1.-]  ) 

D  A  T  A  (  I  r !  X  =  -l,-i,[:,i,i,i,u,-U 

MP  =  J 


.JT  A^iv  ITS  X  A'sl'l 
AKE,  ^LOTTFJ  whcN 


Y  n  I  ^'  f.-  ^1  ■>  lu  N :, 


,  I  in,  IX  0,  I  rn,r.(;'^, 
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CV     .Uf^,     A^C  I  ) 


rv>   ^o   T'^  10 


ISS    =    0 

CV«CL(NLV) 

MTlsMT-i 

Ml    =    MT*l 

DU    10     1=1, Ml 

IF     (AMC I-l)     ,lT. 

IXO     s     I-l 

IX    s     1*1 

lYO     =    1 

lY    =    1 

ISO    =    1 

IS    =    1 

lUX    =    -i 

lOY     =     0 

CALL     TKACt(At1,rx,FY) 
10    COMrjUh 
JsMT-JIil 
DO    20     1=1, Ml 
J    s    J*DIM 

IF     (iM(J«-UIM)     ,lT,     CV     ,ot',     AM(J>     ,nE.     CV)     T  f)     TO    ?0 
IXO     =    MT 
IX    =    .iT 
lYO    =     !♦! 
lY    =    !♦! 
lUX     r     0 
lOY     =     -i 
ISO     =     7 

IS  =  ; 

CALL  TRACt(AM,FA,FY) 
20  COMTNUt 

Jsr^T^NTl*  Jli-1*l 

DO  30  !=l,yTl 

J  =  J-1 

IF  (A.1(J«i)  ,lT,  CV  ,uN,  A'UJ)  ,Gfc,  r;)  nj  Tfi  30 

IXO  s  MT-1 

IX  s  MT-l 

lYO  =  NT 

I  Y  =  NT 

lUX  =  1 

IDY  =  0 

ISO  =  5 

IS  =  t> 

CALL  Tf'Art(AM,rx,FY) 
30  CONTINUb 

JaNT*JI.-l*l 

no  40  1  =  1, Ml 

J  =  J-Dl^i 

IF  ( 4,-i(  j-;jiM)  ,LT,  lv  ,U'<,    am(j)  ,ni,    CV)  '10  to  -iq 


-IS'J- 


ixo    =    1 
IX    =    i 
lYO    =    K'T-i 
lY    =    nT-1 

li>x  =   n 

IDY    =    1 
ISO    =    3 

IS      =      vi 

CALL     TPMC:t(AM,FX,l-'Y) 
40     CUNTIl^lUh 
ISS  =  1 
L    =    0 

no    70    J=2,NT1 
L    =    L^DU" 
DO    oO     1=5  ,Mt 
L    =    L  +  1 

IF     (AH(L*i)     .ur,     CV     ,uK,     Al'lL)     ,nt,     TV)     ^ ']     T'l    ^0 
K    =    L+1 
DQ    bn     iJsl.iJP 
IF    (PtC(  IJ)    .cO,    K)    '.u    lu   ^u 

bO    CONiTI^lLit 
IXO     =     1+3 
IX    =     1*1 
lYO    =    J 
lY    =    J 
IDX    =    -i 

lUY  =  n 

ISO    =    1 

IS    =    i 

CALL    TRACk(AH,I^X,FY) 
60    CfjNTI  Jl't 
70    LsL-MTl 


sijHHnuTir'E    TRAr■l:(A^'.^  A,i- Y) 

UlMb^'SIO^  Arl(l) 

COMhn  >l/l^rLAKC/t^S,l-!0,  ImS,  Tun 

CuMHO  ^l/I^luILE;S/^llH(2)  ,i''l  .NTjKFPL 

CUMHn,N|/XYdNiJS/yA,XM^YA,Yi-!,XbI2F,YbIZP,HX,HY, 
iXS.XSS, YS, YSS,FaA,F YA 

cuMHn,M/CAVIN/DI^  ,  IX,  1  Y,  lux.  inY,  Iss.NP.^^'.rv',  IS,  isn,  IXO,  T  Yn,  PC-='» 
1  I  .^)  X  (  P  )  ,  UJ  Y  (  c  )  ,  f<  u  1  ■  (  ^  u  n  )  ,  X  ( 1  M)  :i )  ,  Y  ( 1  '^  (J 1) 

CUMitPN/GI.  tVtl.b/l  r-L,'!!   v,LL(:3U) 

TrPh     INTF-bFK    hPL.rjIh 
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NxC 

JY  * 

MY  » 

2 

rjaN* 

I^   ( 

!►  ( 

4 

X(N) 

YtN) 

GO  T 

5 

NP  =  w 

PfcC( 

6 

YIN) 

X(fO 

7 

1S=I 

8 

IF  ( 

IS  = 

10 

iux  = 

IDY 

1X2  = 

IY2 

IR  s 

11 

IK  I 

13 

IF(I 

N  s 

X(N) 

Y(N) 

no  T 

15 

IF  ( 

IF  { 

16 

MY  =  D 

IF  ( 

17 

IF  ( 

IX  = 

lY  = 

18 

IS  = 

JY  s 

GO  1 

19 

KYsj 

LY  = 

GO  T 

20 

KY  =  i1 

LY  = 

21 

ncps 

IF  ( 

CALL 

GO  T 

23 

IF  ( 

IX  = 

JIM*(  I  Y-1)*IX 

UlM*IDY*inx*JY 
1 

f.     ,oT,     1600)     GO     To     4U 
IJX)     5,4,6 

=     FlOAT{lY-l)*t-LuAT(I'iY)«(A;i(JY)-CV)/(A'MjY)-Aii(lJlM»irY»JY)) 

s    FLOAK  ly-1 ) 
0    7 
P*l 
'.H)      =     jY 

s    FLOATCIX-I  )*FlOAT{IIU)«(A'-i(JY)-r'V)/(A"(JY)-AM(JY*Iljy)) 

s     I-  L  OAT(  I  Y-1  ) 

IS     ,Lfc,     fl)     GH    TO    lu 

IS-P 
I.>IX(  lb) 
=     I  NY (lb) 

ix*irx 

=     IY*iDY 

IDX^IDY 
S  b  )  1 3  , 1 5 

S,r  t,  ISO.UP.  lY.Nt-.  I  Y0,OH,  IX.NE,  lYQ)     Tj    TO    16 
N*l 

=     X 

s     Y 
0     413 

IX?    ,tO,     0     ,f'R.     lYi^    ,60,     0)    GC    TO    45 
(IXi;     ,GT,     MT)     ,UK,     (lYi?     ,GT,     NT))     GU     T'^     ^5 
IM*IDY>flDX*JY 
IW)     19,17,?0 
CV     ,GT,     AM(r.Y))     (iu    Tl.'    ? 

lY^ 

MY 

n  8 

Y*I  JX 

HY-IUX 
0    21 
Y-IJY 

JY*IJY 
{AM(JY)*Ai^(l\Y)fMlMLY)*A(«1(fY))«,2'' 

cv    ,Lt,    noP)   ijn   ru  2.^ 

GHTPT(AM( JY) ) 

0     7 

\t{     .Gb,     0)     01^     Tu    2^) 
iX2 
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25 


26 


30 
31 

34 

36 

40 
45 


50 
60 

70 
80 

90 

son 


lUX    =     -IPX 

CALL    ljhTPT(  AM(KY)  ) 


=    -ITY 


GtTPK  AIMKY  )  ) 
-IPX 


{  C  V     ,  I  t ,     A  .1  ( 

'  T  (  A  M  (  f-"  Y  )  ) 
^-      n) 


lUY    = 

Gu   m   ?o 
■iy=iy;^ 

lUY     =     -IPY 

CALL 

IX=IXti 

lux   = 

U     (C 
CALL    ufcTP 
IF    (IK    ,nt, 
Ix     =     IX+IUX 
IDX    =    -II'X 
GU    TO    31 
IY=IY*IJY 
lUY    =     -IPY 
U     {  C  V     ,  M  , 
lb    =    IS-1 
JY    =    LY 

nj   TO   lu 

CALL    ^(-TPT(AM(I 
IK     (IK     ,  f  i  t: ,     n  ) 
I  Yt-IUY 
7 


All  (,-iY  )  )     i.ij     TO     10 

Y)) 

bfi    TO    .jP 


h^i  LY  )  )     'ilJ     Tfl    .'.4 


I  Y 

fiU    TO     . 

IX=IY*IUX 

Gu    TO    7 

PRINT    5)Uri, 

II-  (KPUL,F-U 


T(J    Sh. 


U) 


II-  (KPUL,F-U.U)     G 
no    5  0     1=1 ,M 
THfcTA=X( I )*rHb* 
P  =  Y(I)*Kf:-»-RU 


GU     TG     hiJ 


HfcTA  =  X(  I  )*rHb*TNO 

=  Y(  I  )*Kr"»-RU 

(  1  )=XSb*(f^*CUS(THfTA)-|-  XA) 

(  I  )=YSb*(R*SIM(THFTA)-f  YA) 

I'lKll    T  r.ll  ll-; 


P 

X,  .  . 

Y( I )=YS^ 

COM  I -Mil L 
'■•  TO  HO 
DU    .   .     .     _. 

=  XSS*(FX(X  .  .  .  

'  )*YS  + 


GO 

.  _    70 

X  (  I  )      . 

Y(  I  )  =  YSS* 

CUNT INUL 

CaLL    SYMHUL  ( X 

nCJ    9  0     Irl,N 

CALL    PLUT(X( I ) , Y( 


l=i,r' 

"■'    (FX(x(n* 

( F  Y  (  Y  (  I 


X3*XA)-rXA) 
Y  A  )  -  F  Y  A  ) 


,  Y, 


'1 1) 


I' I 


'LV-1,(:,-1  ) 


CONTIrlUfc 
PtTUPN 


\)  ,y  > 


flik 

.     41h 


NAT  (IHu  ,ii3HA     for,  i^u.-     ui. 
■j     IhKMI'lAT(-n     t-'cGAUSt- 


i«ja; 


iT'HjH    LIfh    AT    I.LVtl   ,rl-A 


'>  $ 


IT    LHKTA  iijf-u    -inpr 
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